This file is indexed.

/usr/share/pari/doc/tutorial.tex is in pari-doc 2.5.5-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
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508
2509
2510
2511
2512
2513
2514
2515
2516
2517
2518
2519
2520
2521
2522
2523
2524
2525
2526
2527
2528
2529
2530
2531
2532
2533
2534
2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566
2567
2568
2569
2570
2571
2572
2573
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588
2589
2590
2591
2592
2593
2594
2595
2596
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607
2608
2609
2610
2611
2612
2613
2614
2615
2616
2617
2618
2619
2620
2621
2622
2623
2624
2625
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636
2637
2638
2639
2640
2641
2642
2643
2644
2645
2646
2647
2648
2649
2650
2651
2652
2653
2654
2655
2656
2657
2658
2659
2660
2661
2662
2663
2664
2665
2666
2667
2668
2669
2670
2671
2672
2673
2674
2675
2676
2677
2678
2679
2680
2681
2682
2683
2684
2685
2686
2687
2688
2689
2690
2691
2692
2693
2694
2695
2696
2697
2698
2699
2700
2701
2702
2703
2704
2705
2706
2707
2708
2709
2710
2711
2712
2713
2714
2715
2716
2717
2718
2719
2720
2721
2722
2723
2724
2725
2726
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736
2737
2738
2739
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
2754
2755
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
2776
2777
2778
2779
2780
2781
2782
2783
2784
2785
2786
2787
2788
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799
2800
2801
2802
2803
2804
2805
2806
2807
2808
2809
2810
2811
2812
2813
2814
2815
2816
2817
2818
2819
2820
2821
2822
2823
2824
2825
2826
2827
2828
2829
2830
2831
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841
2842
2843
2844
2845
2846
2847
2848
2849
2850
2851
2852
2853
2854
2855
2856
2857
2858
2859
2860
2861
2862
2863
2864
2865
2866
2867
2868
2869
2870
2871
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882
2883
2884
2885
2886
2887
2888
2889
2890
2891
2892
2893
2894
2895
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929
2930
2931
2932
2933
2934
2935
2936
2937
2938
2939
2940
2941
2942
2943
2944
2945
2946
2947
2948
2949
2950
2951
2952
2953
2954
2955
2956
2957
2958
2959
2960
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
2972
2973
2974
2975
2976
2977
2978
2979
2980
2981
2982
2983
2984
2985
2986
2987
2988
2989
2990
2991
2992
2993
2994
2995
2996
2997
2998
2999
3000
3001
3002
3003
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132
3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
3169
3170
3171
3172
3173
3174
3175
3176
3177
3178
3179
3180
3181
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219
% $Id$
% Copyright (c) 2000  The PARI Group
%
% This file is part of the PARI/\kbd{gp} documentation
%
% Permission is granted to copy, distribute and/or modify this document
% under the terms of the GNU General Public License

% This should be compiled with plain TeX
\def\TITLE{A Tutorial for Pari/GP}
\input parimacro.tex

\chapno=0
\begintitle
\vskip2.5truecm
\centerline{\mine A Tutorial}
\vskip1.truecm
\centerline{\mine for}
\vskip1.truecm
\centerline{\mine PARI / GP}
\vskip1.truecm
\centerline{\sectiontitlebf (version \vers)}
\vskip1.truecm
\authors
\endtitle

\copyrightpage
\tableofcontents

\noindent This booklet is a guided tour and a tutorial to the \kbd{gp}
calculator. Many examples will be given, but each time a new function is
used, the reader should look at the appropriate section in the \emph{User's
Manual to PARI/GP} for detailed explanations. This chapter can be read
independently, for example to get acquainted with the possibilities of
\kbd{gp} without having to read the whole manual. At this point.

\section{Greetings!}

So you are sitting in front of your workstation (or terminal, or PC\dots),
and you type \kbd{gp} to get the program started (or click on the relevant
icon, or select some menu item). It says hello in its particular manner, and
then waits for you after its \kbd{prompt}, initially \kbd{?} (or something
like {\bf gp}~\kbd{>}). Type
\bprog
  2 + 2
@eprog\noindent What happens? Maybe not what you expect. First of all, of
course, you should tell \kbd{gp} that your input is finished, and this is
done by hitting the \kbd{Return} (or \kbd{Newline}, or \kbd{Enter}) key. If
you do exactly this, you will get the expected answer. However some of you
may be used to other systems like Gap, Macsyma, Magma or Maple. In this case,
you will have subconsciously ended the line with a semicolon ``\kbd{;}''
before hitting \kbd{Return}, since this is how it is done on those systems.
In that case, you will simply see \kbd{gp} answering you with a smug
expression, i.e.~a new prompt and no answer!  This is because a semicolon at
the end of a line tells \kbd{gp} not to print the result (it is still stored
in the result history). You will certainly want to use this feature if the
output is several pages long. Try
\bprog
  27 * 37
@eprog\noindent Wow! even multiplication works. Actually, maybe those
spaces are not necessary after all. Let's try \kbd{27*37}. Seems to be ok. We
will still insert them in this document since it makes things easier to read,
but as \kbd{gp} does not care about them, you don't have to type them all.

Now this session is getting lengthy, so the second thing one needs to learn
is to quit. Each system has its quit signal. In \kbd{gp}, you can use
\kbd{quit} or \b{q} (backslash q), the \kbd{q} being of course for quit.
Try it.

Now you've done it! You're out of \kbd{gp}, so how do you want to continue
studying this tutorial? Get back in please.

Let's get to more serious stuff. I seem to remember that the decimal
expansion of $1/7$ has some interesting properties. Let's see what \kbd{gp}
has to say about this. Type
\bprog
  1 / 7
@eprog\noindent
What? This computer is making fun of me, it just spits back to me my own
input, that's not what I want!

Now stop complaining, and think a little. Mathematically, $1/7$ is an element
of the field $\Q$ of rational numbers, so how else but $1/7$ can the computer
give the answer to you? Well maybe $2/14$ or $7^{-1}$, but why complicate
matters? Seriously, the basic point here is that PARI, hence \kbd{gp}, will
almost always try to give you a result which is as precise as possible (we
will see why ``almost'' later). Hence since here the result can be
represented exactly, that's what it gives you.

But I still want the decimal expansion of $1/7$. No problem. Type one of
the following:
\bprog
  1./ 7
  1 / 7.
  1./ 7.
  1 / 7 + 0.
@eprog\noindent
Immediately a number of decimals of this fraction appear, 28 on most systems,
38 on the others, and the repeating pattern is $142857$. The reason is that
you have included in the operations numbers like \kbd{0.}, \kbd{1.} or \kbd{7.}
which are \emph{imprecise} real numbers, hence \kbd{gp} cannot give you an
exact result.

Why 28 / 38 decimals by the way? Well, it is the default initial precision.
This has been chosen so that the computations are very fast, and gives
already 12 decimals more accuracy than conventional double precision floating
point operations. The precise value depends on a technical reason: if your
machine supports 64-bit integers (the standard C library can handle integers
up to $2^{64}$), the default precision is 38 decimals, and 28 otherwise.
For definiteness, we will assume the latter henceforth. Of course, you
can extend the precision (almost) as much as you like as we will see in a
moment.

I'm getting bored, why don't we get on with some more exciting stuff?  Well,
try \kbd{exp(1)}. Presto, comes out the value of $e$ to 28 digits. Try
\kbd{log(exp(1))}. Well, we get a floating point number and not an exact $1$,
but pretty close! That's what you lose by working numerically.

What could we try now? Hum, \kbd{pi}? The answer is not that
enlightening. \kbd{Pi}? Ah. This works better. But let's remember that
\kbd{gp} distinguishes between uppercase and lowercase letters. \kbd{pi} was
as meaningless to it as \kbd{stupid garbage} would have been: in both cases
\kbd{gp} will just create a variable with that funny unknown name you just
used. Try it! Note that it is actually equivalent to type
\kbd{stupidgarbage}: all spaces are suppressed from the input. In the
\kbd{27~*~37} example  it was not so conspicuous as we had an operator to
separate the two operands. This has important consequences for the writing of
\kbd{gp} scripts. More about this later.

By the way, you can ask \kbd{gp} about any identifier you think it might know
about: just type it, prepending a question mark ``\kbd{?}''. Try \kbd{?Pi}
and \kbd{?pi} for instance. On most systems, an extended online help should
be available: try doubling the question mark to check whether it's the case
on yours: \kbd{??Pi}. In fact the \kbd{gp} header already gave you that
information if it was the case, just before the copyright message. As well,
if it says something like ``\kbd{readline enabled}'' then you should have a
look at the \kbd{readline} introduction in the User's Manual before you go
on: it will be much easier to type in examples and correct typos after you've
done that.

Now try \kbd{exp(Pi * sqrt(163))}. Hmmm, we suspect that the last digit may
be wrong, can this really be an integer? This is the time to change
precision. Type \kbd{\b{p} 50}, then try \kbd{exp(Pi * sqrt(163))} again. We
were right to suspect that the last decimal was incorrect, since we get quite
a few nines in its place, but it is now convincingly clear that this is not
an integer. Maybe it's a bug in PARI, and the result is really an integer?
Type
\bprog
  (log(%) / Pi)^2
@eprog\noindent
immediately after the preceding computation; \kbd{\%} means the result of the
last computed expression. More generally, the results are numbered \kbd{\%1,
\%2, \dots} \emph{including} the results
that you do not want to see printed by putting a semicolon at the end of the
line, and you can evidently use all these quantities in any further
computations. The result seems to be indistinguishable from $163$, hence it
does not seem to be a bug.

In fact, it is known that $\exp(\pi*\sqrt{n})$ not only is not an integer or
a rational number, but is even a transcendental number when $n$ is a non-zero
rational number.

So \kbd{gp} is just a fancy calculator, able to give me more decimals than I
will ever need? Not so, \kbd{gp} is incredibly more powerful than an ordinary
calculator, independently of its arbitrary precision possibilities.

\misctitle{Additional comments} (you are supposed to skip this at first,
and come back later)

1) If you are a PARI old timer, say the last version of PARI you used was
released around 1996, you have certainly noticed already that many many
things changed between the older 1.39.xx versions and this one.
Conspicuously, most function names have been changed.

Of course, this is going to break all your nice old scripts. Well, you can
either change the compatibility level (typing \kbd{default(compatible, 3)}
will send you back to the stone-age behavior of good ol' version 1.39.15),
or rewrite the scripts. We really advise you to do the latter if they are not
too long, since they can now be written much more cleanly than before
(especially with the new control statements: \kbd{break}, \kbd{next},
\kbd{return}), and besides it'll be as good a way as any to get used to the
new names.

To know how a specific function was changed, just type \kbd{whatnow({\rm
function})}.

2) It seems that the text implicitly says that as soon as an imprecise number
is entered, the result will be imprecise. Is this always true? There is a
unique exception: when you multiply an imprecise number by the exact number
0, you will get the exact 0. Compare \kbd{0 * 1.4} and \kbd{0.~*~1.4}.
\smallskip
%
3) Not only can the number of decimal places of real numbers be large, but
the number of digits of integers also. Try \kbd{1000!}. It is never necessary
to tell \kbd{gp} in advance the size of the integers that it will encounter.
The same is true for real numbers, although most computations with floating
point assume a default precision and truncate their results to this accuracy;
initially 28 decimal digits, but we may change that with \b{p} of course.
\smallskip
%
4) Come back to 28 digits of precision (\kbd{\b{p} 28}), and type
\kbd{exp(75)}. As you can see the result is printed in exponential format.
This is because \kbd{gp} never wants you to believe that a result is correct
when it is not. We are working with 28 digits of precision, but the integer
part of $\exp(24*\pi)$ has 33 decimal digits. Hence if \kbd{gp} had dutifully
printed out 33 digits, the last few digits would have been wrong. Hence
\kbd{gp} wants to print only 28 significant digits, but to do so it has to
print in exponential format. \smallskip
%
5) There are two ways to avoid this. One is of course to increase the
precision to more than 33 decimals. Let's try it. To give it a wide margin,
we set the precision to 40 decimals. Then we recall our last result (\kbd{\%}
or \kbd{\%n} where \kbd{n} is the number of the result). What? We still have
an exponential format! Do you understand why?

Again let's try to see what's happening. The number you recalled had been
computed only to 28 decimals, and even if you set the precision to 1000
decimals, \kbd{gp} knows that your number has only 28 digits of accuracy but
an integral part with 33 digits. So you haven't improved things by increasing
the precision. Or have you? What if we retype \kbd{exp(75)} now that we
have 40 digits? Try it. Now we do not have an exponential format, and we see
that at 28 decimals the last 6 digits would have been wrong if they had been
printed in fixed-point format.
\medskip
%
6) What if I forget what the current precision is and I don't feel like
counting all the decimals? Well, you can type \b{p} by itself. You may also
learn about \kbd{gp} internal variables (and change them!) using
\kbd{default}. Type \kbd{default(realprecision)}, then
\kbd{default(realprecision, 38)}. Huh? In fact this last command is strictly
equivalent to \kbd{\b{p} 38}! (Admittedly more cumbersome to type.) There are
more ``defaults'' than just \kbd{format} and \kbd{realprecision}: type
\kbd{default} by itself now, they are all there. \smallskip
%
7) Note that the \kbd{default} command reacts differently according to the
number of input arguments. This is not an uncommon behavior for \kbd{gp}
functions. You can see this from the online help, or the complete description
in Chapter~3: any argument surrounded by braces \kbd{\obr\cbr} in the
function prototype is optional, which really means that a \emph{default}
argument will be supplied by \kbd{gp}. You can then check out from the text
what effect a given value will have, and in particular the default one.
\smallskip
%
8) Try the following: starting in precision 28, type first
\kbd{default(format, "e0.100")}, then \kbd{exp(1)}. Where are my 100
significant digits? Well, \kbd{default(format,)} only changes the output
format, but \emph{not} the default precision. On the other hand, the \b{p}
command changes both the precision and the output format.

\section{Warming up}

Another thing you better get used to pretty fast is error messages. Try
typing \kbd{1/0}. Couldn't be clearer. Taking again our universal example in
precision 28, type
\bprog
  floor(exp(75))
@eprog\noindent
\kbd{floor} is the mathematician's integer part, not to be confused with
\kbd{truncate}, which is the computer scientist's: \kbd{floor(-3.4)} is equal
to $-4$ whereas \kbd{truncate(-3.4)} is equal to $-3$.  You get a more
cryptic error message, which you would immediately understand if you had read
the additional comments of the preceding section. Since you were told not to
read them, here's the explanation: \kbd{gp} is unable to compute the
integer part of \kbd{exp(75)} given only 28 decimals of accuracy, since
it has 33 digits.

Some error messages are more cryptic and sometimes not so easy to understand.
For instance, try \kbd{log(x)}. It simply tells you that \kbd{gp} does not
understand what \kbd{log(x)} is, although it does know the \kbd{log}
function, as \kbd{?log} will readily tell us.

Now let's try \kbd{sqrt(-1)} to see what error message we get now. Haha!
\kbd{gp} even knows about complex numbers, so impossible to trick it that
way. Similarly, try typing \kbd{log(-2)}, \kbd{exp(I*Pi)}, \kbd{I\pow
I}\dots\ So we have a lot of real and complex analysis at our disposal.
There always is a specific branch of multivalued complex transcendental
functions which is taken, specified in the manual. Again, beware that
\kbd{I} and \kbd{i} are not the same thing. Compare \kbd{I\pow2} with
\kbd{i\pow2} for instance.

Just for fun, let's try \kbd{6*zeta(2) / Pi\pow2}. Pretty close, no?

\medskip
Now \kbd{gp} didn't seem to know what \kbd{log(x)} was, although it did know
how to compute numerical values of \kbd{log}. This is annoying. Maybe it
knows the exponential function? Let's give it a try. Type \kbd{exp(x)}.
What's this? If you had any experience with other computer algebra systems,
the answer should have simply been \kbd{exp(x)} again. But here the answer is
the Taylor expansion of the function around $\kbd{x}=0$, to 16 terms. 16 is
the default \kbd{seriesprecision}, which can be changed by typing \kbd{\b{ps}
$n$} or \kbd{default(seriesprecision, $n$)} where $n$ is the number of terms
that you want in your power series. Note the \kbd{O(x\pow16)} which ends the
series, and which is trademark of this type of object in \kbd{gp}. It is the
familiar ``big--oh'' notation of analysis.

You thus automatically get the Taylor expansion of any function that can be
expanded around $0$, and incidentally this explains why we weren't able to do
anything with \kbd{log(x)} which is not defined at $0$. (In fact \kbd{gp}
knows about Laurent series, but \kbd{log(x)} is not meromorphic either at
$0$.) If we try \kbd{log(1+x)}, then it works. But what if we wanted the
expansion around a point different from 0? Well, you're able to change $x$
into $x-a$, aren't you? So for instance you can type \kbd{log(x+2)} to have
the expansion of \kbd{log} around $\kbd{x}=2$. As exercises you can try
\bprog
  cos(x)
  cos(x)^2 + sin(x)^2
  exp(cos(x))
  gamma(1 + x)
  exp(exp(x) - 1)
  1 / tan(x)
@eprog\noindent
for different values of \kbd{serieslength} (change it using \b{ps}
\var{newvalue}).

Let's try something else: type \kbd{(1 + x)\pow 3}. No \kbd{O(x)} here, since
the result is a polynomial.  Haha, but I have learnt that if you do not take
exponents which are integers greater or equal to 0, you obtain a power series
with an infinite number of non-zero terms. Let's try.  Type
\kbd{(1 + x)\pow (-3)} (the parentheses around \kbd{-3} are not necessary but
make things easier to read). Surprise! Contrary to what we expected, we don't
get a power series but a rational function. Again this is for the same reason
that \kbd{1 / 7} just gave you $1/7$: the result being exact, PARI doesn't see
any reason to make it non-exact.

But I still want that power series. To obtain it, you can do as in the $1/7$
example and type
\bprog
  (1 + x)^(-3) + O(x^16)
  (1 + x)^(-3) * (1 + O(x^16))
  (1 + x + O(x^16))^(-3)
@eprog\noindent
(Not on this example, but there is a difference between the first $2$
methods. Do you spot it?) Better yet, use the series constructor which
transforms any object into a power series, using the current
\kbd{seriesprecision}, and simply type
\bprog
  Ser( (1 + x)^(-3) )
@eprog

Now try \kbd{(1 + x)\pow (1/2)}: we obtain a power series, since the
result is an object which PARI does not know how to represent exactly. (We
could teach PARI about algebraic functions, but then take \kbd{(1 + x)\pow Pi}
as another example.) This gives us still another solution to our preceding
exercise: we can type \kbd{(1 + x)\pow (-3.)}. Since \kbd{-3.} is not an exact
quantity, PARI has no means to know that we are dealing with a rational
function, and will instead give you the power series, this time with real
instead of integer coefficients.
\smallskip

To summarize, in this section we have seen that in addition to integers, real
numbers and rational numbers, PARI can handle complex numbers, polynomials,
rational functions and power series. A large number of functions exist which
handle these types, but in this tutorial we will only look at a few.

\misctitle{Additional comments} (as before, you are supposed to skip this
at first reading)

1) A complex number has a real part and an imaginary part (who would have
guessed?). However, beware that when the imaginary part is the exact integer
zero, it is not printed, but the complex number is not converted to a real
number. Hence it may \emph{look} like a real number without being one, and
this may cause some confusion in programs which expect real numbers. For
example, type \kbd{n = 3 + 0*I}. The answer reads \kbd{3}, as expected. But
it is still a complex number. Check it with \kbd{type(n)}. Hence if you then
type \kbd{(1+x)\pow n}, instead of getting the polynomial \kbd{1 + 3*x +
3*x\pow 2 + x\pow 3} as expected, you obtain a power series. Worse, when you
try to apply an arithmetic function, say the Euler totient function (known as
\kbd{eulerphi} to \kbd{gp}), you get an error message worrying about integer
arguments. You would have guessed yourself, but the message is difficult to
understand since 3 looks like a genuine integer! Please read again if this
is not clear; it is a common source of incomprehension.

Similarly, \kbd{n = 3 + 0*x} is not the integer 3 but a constant polynomial
equal to $3x^0$.

2) If you want the final expression to be in the simplest form possible (for
example before applying an arithmetic function, or simply because things will
go faster afterwards), apply the function \kbd{simplify} to the result.
This is done automatically at the end of a \kbd{gp} command, but
\emph{not} in intermediate expressions. Hence \kbd{n} above is not an
integer, but the final result stored in the output history is! So
if you type \kbd{type(\%)} instead of \kbd{type(n)} the answer is
\typ{INT}, adding to the confusion.

3) As already stated, power series expansions are always implicitly around
$\kbd{x} = 0$. When we wanted them around $\kbd{x} = \kbd{a}$, we replaced
\kbd{x} by \kbd{z + a} in the function we wanted to expand. For complicated
functions, it may be simpler to use the substitution function \kbd{subst}.
For example, if \kbd{p~= 1 / (x\pow 4 + 3*x\pow 3 + 5*x\pow 2 - 6*x + 7)},
you may not want to retype this, replacing \kbd{x} by \kbd{z~+ a}, so you can
write \kbd{subst(p, x, z+a)} (look up the exact description of the
\kbd{subst} function).

Now type \kbd{subst(1 + O(x), x, z+1)}. Do you understand the error message?

4) The valuation at $\kbd{x} = 0$ for a power series \kbd{p} is obtained
as \kbd{valuation(p, x)}.

\section{The Remaining PARI Types}
Let's talk some more about the basic PARI types.

Type \kbd{p = x * exp(-x)}. As expected, you get the power series expansion
to 16 terms (if you have not changed the default). Now type
\kbd{pr = serreverse(p)}. You are asking here for the \emph{reversion} of the
power series \kbd{p}, in other words the inverse function. This is possible
only for power series whose first non-zero coefficient is that of $x^1$.  To
check the correctness of the result, you can type \kbd{subst(p, x, pr)} or
\kbd{ subst(pr, x, p)} and you should get back \kbd{x + O(x\pow 17)}.

Now the coefficients of \kbd{pr} obey a very simple formula. First, we would
like to multiply the coefficient of \kbd{x\pow n} by \kbd{n!} (in the case of
the exponential function, this would simplify things considerably!). The PARI
function \kbd{serlaplace} does just that. So type \kbd{ps = serlaplace(pr)}.
The coefficients now become integers, which can be immediately recognized by
inspection. The coefficient of $x^n$ is now equal to
$n^{n-1}$. In other words, we have
%
$$\kbd{pr} = \sum_{n\ge1}\dfrac{n^{n-1}}{n!} X^{n}.$$
%
Do you know how to prove this? (The proof is difficult.)
\smallskip
%
Of course PARI knows about vectors (rows and columns are distinguished, even
though mathematically there is no difference) and matrices. Type for example
\kbd{[1,2,3,4]}. This gives the row vector whose coordinates are 1, 2, 3 and
4.  If you want a column vector, type \kbd{[1,2,3,4]\til}, the tilde meaning
of course transpose. You don't see much difference in the output, except for
the tilde at the end. However, now type \b{b}: lo and behold, the column
vector appears as a proper vertical thingy now. The \b{b} command is used
mainly for this purpose. The length of a vector is given by, well
\kbd{length} of course. The shorthand ``cardinal'' notation \kbd{\#v} for
\kbd{length(v)} is also available, for instance \kbd{v[\#v]} is the last
element of \kbd{v}.

Type \kbd{m = [a,b,c; d,e,f]}. You have just entered a matrix with 2 rows and
3 columns. Note that the matrix is entered by \emph{rows} and the rows are
separated by semicolons ``\kbd{;}''. The matrix is printed naturally in a
rectangle shape. If you want it printed horizontally just as you typed it,
type \b{a}, or if you want this type of printing to be the permanent default
type \kbd{default(output, 0)}. Type \kbd{default(output, 1)} if you want to
come back to the original output mode.

Now type \kbd{m[1,2]}, \kbd{m[1,]}, \kbd{m[,2]}. Are explanations necessary?
(In an expression such as \kbd{m[j,k]}, the \kbd{j} always refers to the
row number, and the \kbd{k} to the column number, and the first index is
always 1, never 0. This default cannot be changed.)

Even better, type \kbd{m[1,2] = 5; m}. The semicolon also allows us to put
several instructions on the same line; the final result is the output of
the last statement on the line. Now type \kbd{m[1,] = [15,-17,8]}. No
problem. Finally type \kbd{m[,2] = [j,k]}. You have an error message since you
have typed a row vector, while \kbd{m[,2]} is a column vector. If you type
instead \kbd{m[,2] = [j,k]\til} it works. \smallskip
%
\label{se:types}
Type now \kbd{h = mathilbert(20)}. You get the so-called ``Hilbert matrix''
whose coefficient of row $i$ and column $j$ is equal to $(i+j-1)^{-1}$.
Incidentally, the matrix \kbd{h} takes too much room. If you don't want to
see it, simply type a semi-colon ``\kbd{;}'' at the end of the line
(\kbd{h = mathilbert(20);}). This is an example of a ``precomputed'' matrix,
built into PARI. We will see a more general construction later.

What is interesting about Hilbert matrices is that first their inverses and
determinants can be computed explicitly (and the inverse has integer
coefficients), and second they are numerically very unstable, which make them
a severe test for linear algebra packages in numerical analysis.  Of course
with PARI, no such problem can occur: since the coefficients are given as
rational numbers, the computation will be done exactly, so there cannot be
any numerical error. Try it. Type \kbd{d~=~matdet(h)}. The result is a
rational number (of course) of numerator equal to 1 and denominator having
226 digits. How do I know, by the way? Well, type \kbd{sizedigit(1/d)}. Or
\kbd{\#Str(1/d)}. (The length of the character string representing the
result.)

Now type \kbd{hr = 1.* h;} (do not forget the semicolon, we don't want to see
the result!), then \kbd{dr = matdet(hr)}. You notice two things. First the
computation, is much faster than in the rational case. (If your computer is
too fast for you to notice, try again with \kbd{h = mathilbert(40)}, or
some larger value.) The reason for this is that PARI is handling real
numbers with 28 digits of accuracy, while in the rational case it is
handling integers having up to 226 decimal digits.

The second, more important, fact is that the result is terribly wrong. If you
compare with \kbd{1.$*$d} computed earlier, which is the correct answer, you
will see that only 2 decimals agree! (None agree if you replaced 20 by 40 as
suggested above.) This catastrophic instability is as already mentioned one
of the characteristics of Hilbert matrices. In fact, the situation is
worse than that. Type \kbd{norml2(1/h - 1/hr)} (the function \kbd{norml2}
gives the square of the $L^2$ norm, i.e.~the sum of the squares of the
coefficients). The result is larger than $10^{50}$, showing that some
coefficients of \kbd{1/hr} are wrong by as much as $10^{24}$ (the largest
error is in fact equal to $4.244 \cdot 10^{24}$ for the coefficient of row 15
and column 15, which is a 28 digit integer). To obtain the correct result
after rounding for the inverse, we have to use a default precision of 56
digits (try it).

Although vectors and matrices can be entered manually, by typing explicitly
their elements, very often the elements satisfy a simple law and one uses a
different syntax. For example, assume that you want a vector whose $i$-th
coordinate is equal to $i^2$. No problem, type for example
\kbd{vector(10,i, i\pow 2)} if you want a vector of length 10. Similarly, if
you type
\bprog
  matrix(5,5, i,j, 1 / (i+j-1))
@eprog\noindent
you will get the Hilbert matrix of order 5, hence the \kbd{mathilbert}
function is in fact redundant.  The \kbd{i} and \kbd{j} represent dummy
variables which are used to number the rows and columns respectively (in
the case of a vector only one is present of course). You must not forget,
in addition to the dimensions of the vector or matrix, to indicate
explicitly the names of these variables. You may omit the variables and
the final expression to get zero entries, as in \kbd{matrix(10,20)}.

\misctitle{Warning} The letter \kbd{I} is reserved for the complex number
equal to the square root of $-1$. Hence it is forbidden to use it as a
variable. Try typing \kbd{vector(10,I, I\pow 2)}, the error message that you
get clearly indicates that \kbd{gp} does not consider \kbd{I} as a variable.
There are two other reserved variable names: \kbd{Pi} and \kbd{Euler}. All
function names are forbidden as well. On the other hand there is nothing
special about \kbd{i}, \kbd{pi} or \kbd{euler}.

When creating vectors or matrices, it is often useful to use Boolean
operators and the \kbd{if()} statement. Indeed, an \kbd{if} expression has a
value, which is of course equal to the evaluated part of the \kbd{if}. So for
example you can type
\bprog
  matrix(8,8, i,j, if ((i-j)%2, 1, 0))
@eprog\noindent
to get a checkerboard matrix of \kbd{1} and \kbd{0}. Note however
that a vector or matrix must be \emph{created} first before being used. For
example, it is possible to write
\bprog
  v = vector(5);
  for (i = 1, 5, v[i] = 1/i)
@eprog\noindent
but this would fail if the vector \kbd{v} had not been created beforehand.
Another useful way to create vectors and matrices is to extract them from
larger ones, using \kbd{vecextract}. For instance, if \kbd{h} is the
$20\times 20$ Hilbert matrix as above,
\bprog
  h = mathilbert(20);
  vecextract(h, "11..20", "11..20")
@eprog\noindent is its lower right quadrant.

\medskip The last PARI types which we have not yet played with are closely
linked to number theory. People not interested in number theory can skip
ahead.

The first is the type ``integer--modulo''. Let us see an example. Type
\bprog
  n = 10^15 + 3
@eprog
We want to know whether this number is prime or not. Of course we could make
use of the built-in facilities of PARI, but let us do otherwise. We first
trial divide by the built-in table of primes. We slightly cheat here and use
a variant of the function \kbd{factor} which does exactly this. So type
\kbd{factor(n, 200000)}. The last argument tells \kbd{factor} to trial divide
up to the given bound and stop at this point. Set it to 0 to trial divide by
the full set of built-in primes, which goes up to $500000$ by default.

As for all factoring functions, the result is a 2 column matrix: the first
column gives the primes and the second their exponents. Here we get a single
row, telling us that if primes stopped at $200000$ as we made \kbd{factor}
believe, \kbd{n} would be prime. (Or is that a contradiction?) More
seriously, \kbd{n} is not divisible by any prime up to $200000$.

We could now trial divide further, or cheat and call the PARI function
\kbd{factor} without the optional second argument, but before we do this let
us see how to get an answer ourselves.

By Fermat's little theorem, if $n$ is prime we must have $a^{n-1}\equiv 1
\pmod{n}$ for all $a$ not divisible by $n$. Hence we could try this with $a=2$
for example. But $2^{n-1}$ is a number with approximately $3\cdot10^{14}$
digits, hence impossible to write down, let alone to compute. But instead type
\kbd{a = Mod(2,n)}. This creates the number $2$ considered now as an element
of the ring $R = \Z/\kbd{n}\Z$. The elements of $R$, called intmods, can
always be represented by numbers smaller than \kbd{n}, hence small. Fermat's
theorem can be rewritten
%
$\kbd{a}^{n-1} = \kbd{Mod(1,n)}$
%
in the ring $R$, and this can be computed very efficiently. Elements of $R$
may be lifted back to $\Z$ with either \kbd{lift} or \kbd{centerlift}. Type
\kbd{a\pow (n-1)}. The result is definitely \emph{not} equal to
\kbd{Mod(1,n)}, thus \emph{proving} that \kbd{n} is not a prime. If we had
obtained \kbd{Mod(1,n)} on the other hand, it would have given us a hint that
\kbd{n} is maybe prime, but never a proof.

To find the factors is another story. One must use less naive techniques than
trial division, or be patient. To finish this example, type \kbd{fa =
factor(n)} to see the factors. The smallest factor is 14902357, we may take
a small coffee break and let trial division run to completion. Type
\kbd{\#} to turn the time on, then
\bprog
  for (i = 2, ceil(sqrt(n)), if (n%i==0, print(i); break))
@eprog\noindent
This takes less than one minute on a 1GHz processor. Actually, interpreting
the expression in the loop takes more than two thirds of the running time
as you can check by typing just
\bprog
  for (i = 2, 14902357, n%i)
@eprog\noindent
You may stop the timer by typing \kbd{\#} again.

Note that, as is the case with most ``prime''-producing functions, the
``prime'' factors given by \kbd{factor} are only strong pseudoprimes, and not
\emph{proven} primes.  Use \kbd{isprime( fa[,1] )} to rigorously prove
primality of the factors. The latter command applies \kbd{isprime} to all
entries in the first column of \kbd{fa}, i.e to all pseudoprimes, and returns
the column vector of results: all equal to 1, so our pseudoprimes were
true primes. All arithmetic functions can be applied in this way to the entries
of a vector or matrix. In fact, it has been checked that the strong
pseudoprimes output by \kbd{factor} (Baillie-Pomerance-Selfridge-Wagstaff
pseudoprimes, without small divisors) are true primes at least up to
$10^{13}$, and no explicit counter-example is known.\smallskip

The second specifically number-theoretic type is the $p$-adic numbers. I have
no room for definitions, so please skip ahead if you have no use for such
beasts. A $p$-adic number is entered as a rational or integer valued
expression to which is added \kbd{O(p\pow n)}, or simply \kbd{O(p)} if
$\kbd{n}=1$, where \kbd{p} is the prime and \kbd{n} the $p$-adic precision.
Note that you have to explicitly type in \kbd{3\pow 2} for instance, \kbd{9}
will not do. Unless you want to cheat \kbd{gp} into believing that \kbd{9}
is prime, but you had better know what you are doing in this case: most
computations will yield a wrong result.

Apart from the usual arithmetic operations, you can apply a number of
transcendental functions. For example, type \kbd{n = 569 + O(7\pow 8)}, then
\kbd{s~=~sqrt(n)}, you obtain one of the square roots of \kbd{n}; to check
this, type \kbd{s\pow 2 - n}). Type now \kbd{s = log(n)}, then \kbd{e =
exp(s)}. If you know about $p$-adic logarithms, you will not be surprised
that \kbd{e} is not equal to \kbd{n}. Type \kbd{(n/e)\pow 6}: \kbd{e} is in
fact equal to \kbd{n} times the $(p-1)$-st root of unity \kbd{teichmuller(n)}.

Incidentally, if you want to get back the integer 569 from the $p$-adic
number \kbd{n}, type \kbd{lift(n)} or \kbd{truncate(n)}.
\smallskip

The third number-theoretic type is the type ``quadratic number''. This type
is specially tailored so that we can easily work in a quadratic extension of
a base field, usually $\Q$. It is a generalization of the type
``complex''. To start, we must specify which quadratic field we want to work
in. For this, we use the function \kbd{quadgen} applied to the
\emph{discriminant} \kbd{d} (as opposed to the radicand) of the quadratic
field. This returns a number (always printed as \kbd{w}) equal to
$(\kbd{d}+a) / 2$ where $a$ is equal to 0 or 1 according to whether \kbd{d} is
even or odd. The behavior of \kbd{quadgen} is a little special: although its
result is always printed as \kbd{w}, the variable \kbd{w} itself is not set
to that value. Hence it is necessary to write systematically
\kbd{w = quadgen(d)} using the variable name \kbd{w} (or \kbd{w1} etc. if you
have several quadratic fields), otherwise things will get confusing.

So type \kbd{w = quadgen(-163)}, then \kbd{charpoly(w)} which asks for the
characteristic polynomial of \kbd{w}. The result shows what \kbd{w} will
represent. You may ask for \kbd{1.*w} to see which root of the quadratic has
been taken, but this is rarely necessary. We can now play in the field
$\Q(\sqrt{-163})$. Type for example \kbd{w\pow 10}, \kbd{norm(3 + 4*w)},
\kbd{1 / (4+w)}. More interesting, type \kbd{a = Mod(1,23) * w} then \kbd{b =
a\pow 264}. This is a generalization of Fermat's theorem to quadratic fields.
If you do not want to see the modulus 23 all the time, type \kbd{lift(b)}.

Another example: type \kbd{p = x\pow 2 + w*x + 5*w + 7}, then \kbd{norm(p)}. We
thus obtain the quartic equation over $\Q$ corresponding to the relative
quadratic extension over $\Q(\kbd{w})$ defined by \kbd{p}.

On the other hand, if you type \kbd{wr  = sqrt(w\pow 2)}, do not expect to get
back \kbd{w}. Instead, you get the numerical value, the function \kbd{sqrt}
being considered as a ``transcendental'' function, even though it is
algebraic. Type \kbd{algdep(wr,2)}: this looks for algebraic relations
involving the powers of \kbd{w} up to degree 2. This is one way to get
\kbd{w} back. Similarly, type \kbd{algdep(sqrt(3*w + 5), 4)}. See the user's
manual for the function \kbd{algdep}.\smallskip

The fourth number-theoretic type is the type ``polynomial--modulo'', i.e.
polynomial modulo another polynomial. This type is used to work in general
algebraic extensions, for example elements of number fields (if the base
field is $\Q$), or elements of finite fields (if the base field is
$\Z/p\Z$ for a prime $p$). In a sense it is a generalization of the type
quadratic number. The syntax used is the same as for intmods. For example,
instead of typing \kbd{w = quadgen(-163)}, you can type
\bprog
  w = Mod(x, quadpoly(-163))
@eprog\noindent
Then, exactly as in the quadratic case, you can type \kbd{w\pow 10},
\kbd{norm(3 + 4*w)}, \kbd{1 / (4+w)}, \kbd{a = Mod(1,23)*w}, \kbd{b = a\pow
264}, obtaining of course the same results. (Type \kbd{lift(\dots)} if you
don't want to see the polynomial \kbd{x\pow 2 - x + 41} repeated all the
time.) Of course, you can work in any degree, not only quadratic. For the
latter, the corresponding elementary operations will be slower than
with quadratic numbers. Start the timer, then compare
\bprog
  w = quadgen(-163); W = Mod(x, quadpoly(-163));
  a = 2 + w;         A = 2 + W;
  b = 3 + w;         B = 3 + W;
  for (i=1,10^5, a+b)
  for (i=1,10^5, A+B)
  for (i=1,10^5, a*b)
  for (i=1,10^5, A*B)
  for (i=1,10^5, a/b)
  for (i=1,10^5, A/B)
@eprog\noindent
Don't retype everything, use the arrow keys!

There is however a slight difference in behavior. Keeping our polmod \kbd{w},
type \kbd{1.*w}. As you can see, the result is not the same. Type
\kbd{sqrt(w)}. Here, we obtain a vector with 2 components, the two components
being the principal branch of the square root of all the possible embeddings
of \kbd{w} in $\C$. More generally, if
\kbd{w} was of degree $n$, we would get an $n$-component vector, and similarly
for all transcendental functions.

We have at our disposal the usual arithmetic functions, plus a few others.
Type \kbd{a = Mod(x, x\pow 3 - x - 1)} defining a cubic extension. We can for
example ask for \kbd{b = a\pow 5}. Now assume we want to express \kbd{a}
as a polynomial in \kbd{b}. This is possible since \kbd{b} is also a
generator of the same field. No problem, type \kbd{modreverse(b)}. This gives
a new defining polynomial for the same field, i.e.~the characteristic
polynomial of \kbd{b}, and expresses \kbd{a} in terms of this new polmod,
i.e.~in terms of \kbd{a}. We will see this in more detail in the number
field section.

\section{Elementary Arithmetic Functions}

Since PARI is aimed at number theorists, it is not surprising that there
exists a large number of arithmetic functions; see the list by typing
\kbd{?4}. We have already seen several, such as \kbd{factor}. Note that
\kbd{factor} handles not only integers, but also univariate polynomials.
Type for example \kbd{factor(x\pow 200 - 1)}. You can also ask to factor a
polynomial modulo a prime $p$ (\kbd{factormod}) and even in a finite field
which is not a prime field (\kbd{factorff}).

Evidently, you have functions for computing GCD's (\kbd{gcd}), extended GCD's
(\kbd{bezout}), solving the Chinese remainder theorem (\kbd{chinese}) and so
on.

In addition to the factoring facilities, you have a few functions related to
primality testing such as \kbd{isprime}, \kbd{ispseudoprime},
\kbd{precprime}, and \kbd{nextprime}. As previously mentioned, only strong
pseudoprimes are produced by the latter two (they pass the
\kbd{ispseudoprime} test); the more sophisticated primality tests in
\kbd{isprime}, being so much slower, are not applied by default.

We also have the usual multiplicative arithmetic functions: the M\"obius $\mu$
function (\kbd{moebius}), the Euler $\phi$ function (\kbd{eulerphi}), the
$\omega$ and $\Omega$ functions (\kbd{omega} and \kbd{bigomega}), the
$\sigma_k$ functions (\kbd{sigma}), which compute sums of $k$-th powers of the
positive divisors of a given integer, etc\dots

You can compute continued fractions. For example, type \kbd{\b{p} 1000}, then
\kbd{contfrac(exp(1))}: you obtain the continued fraction of the base of
natural logarithms, which as you can see obeys a very simple pattern. Can
you prove it?

In many cases, one wants to perform some task only when an arithmetic
condition is satisfied. \kbd{gp} gives you the following functions: \kbd{isprime}
as mentioned above, \kbd{Z\_issquare}, \kbd{isfundamental} to test whether an
integer is a fundamental discriminant (i.e.~$1$ or the discriminant of a
quadratic field), and the \kbd{forprime}, \kbd{fordiv} and \kbd{sumdiv}
loops. Assume for example that we want to compute the product of all the
divisors of a positive integer \kbd{n}. The easiest way is to write
\bprog
  p = 1; fordiv(n,d, p *= d); p
@eprog\noindent
(There is a simple formula for this product in terms of $n$ and the number of
its divisors: find and prove it!) The notation \kbd{p *= d} is just a
shorthand for \kbd{p = p * d}.

If we want to know the list of primes $p$ less than 1000 such that 2 is a
primitive root modulo $p$, one way would be to write:
\bprog
  forprime(p=3,1000, if (znprimroot(p) == 2, print(p)))
@eprog\noindent
%
Note that this assumes that \kbd{znprimroot} returns the smallest primitive
root, and this is indeed the case. Had we not known about this, we could
have written
\bprog
  forprime(p=3,1000, if (znorder(Mod(2,p)) == p-1, print(p)))
@eprog\noindent
%
(which is actually faster since we only compute the order of $2$ in $\Z/p\Z$,
instead of looking for a generator by trying successive elements whose orders
have to be computed as well.)

Arithmetic functions related to quadratic fields, binary quadratic forms and
general number fields will be seen in the next sections.

\section{Performing Linear Algebra}
The standard linear algebra routines are available: \kbd{matdet},
\kbd{mateigen} (eigenvectors), \kbd{matker}, \kbd{matimage}, \kbd{matrank},
\kbd{matsolve} (to solve a linear system), \kbd{charpoly} (characteristic
polynomial), to name a few. Bilinear algebra over $\R$ is also there:
\kbd{qfgaussred} (Gauss reduction), \kbd{qfsign} (signature). You may also
type \kbd{?8}. Can you guess what each of these do?

Let us see how this works. First, a vector space (or module) is given by a
generating set of vectors (often a basis) which are represented as
\emph{column} vectors. This set of vectors is in turn represented by the
columns of a matrix. Quadratic forms are represented by their Gram matrix.
The base field (or ring) can be any ring type PARI supports. However, certain
operations are specifically written for a real or complex base field, while
others are written for $\Z$ as the base ring.

We had some fun with Hilbert matrices and numerical instability a while back,
but most of the linear algebra routines are generic. If as before \kbd{h =
mathilbert(20)}, we may compute
\bprog
  matdet(h * Mod(1,101))
  matdet(h * (1 + O(101^100)))
@eprog\noindent
in $\Z/101\Z$ and the $p$-adic ring $\Z_{101}$ (to $100$ words of accuracy)
respectively. Let \kbd{H = 1/h} the inverse of \kbd{h}:
\bprog
  H = 1/h;  \\ @com integral
  L = vecextract( primes(10000), "-50.." );  \\ @com extract the last 50 elements
  v = vector(#L, i, matdet(H * Mod(1,L[i])));
  centerlift( chinese(v) )
@eprog\noindent
returns the determinant of \kbd{H}. (Assuming it is an integer
less than half the product of elements of \kbd{L} in absolute value, which
it is.)
In fact, we computed an homomorphic image of the determinant in a few small
finite fields, which admits a single integer representative given the size
constraints. We could also have made a single determinant computation modulo
a big prime (or pseudoprime) number, e.g \kbd{nextprime(2 * B)} if we know
that the determinant is less than \kbd{B} in absolute value.
(Why is that $2$ necessary?)

By the way, this is how you insert comments in a script: everything
following a double backslash, up to the first newline character, is ignored.
If you want comments which span many lines, you can brace them between
\kbd{/* ... */} pairs. Everything in between will be ignored as well. For
instance as a header for the script above you could insert the
following:
\bprog
  /* Homomorphic imaging scheme to compute the determinant of a classical
   * integral matrix.
   * TODO: Look up the explicit formula
   */
@eprog\noindent
(I hope you did not waste your time copying this nonsense, did you?)
\medskip

In addition, linear algebra over $\Z$, i.e.~work on lattices, can also be
performed. Let us now consider the lattice $\Lambda$ generated by the columns
of \kbd{H} in $\Z^{20}\subset\R^{20}$. Since the determinant is non-zero, we
have in fact a basis. What is the structure of the finite abelian group
$\Z^{20}/\Lambda$? Type \kbd{matsnf(H)}. Wow, 20 cyclic factors.


There is a triangular basis for $\Lambda$ (triangular when expressed in
the canonical basis), perhaps it looks better than our initial one? Type
\kbd{mathnf(H)}. Hum. This is taking ages. Ah, the documentation says that
\kbd{mathnf} is fast for tiny matrices but should never be used for matrices
of dimension larger than $10$. Type \kbd{Control-C} to interrupt the longish
computation, then \kbd{z = mathnf(H, 1);}. That's better. \kbd{z[1]} is the
triangular HNF basis, and \kbd{z[2]} is the base change matrix from the
canonical basis to the new one; of course it has determinant $\pm 1$. Try
\kbd{matdet(z[2])}, then \kbd{H * z[2] == z[1]}. Fine, it works. And
\kbd{z[1]} indeed looks better than \kbd{H}.

Can we do better? Perhaps, but then we'd better drop the requirement that
the basis be triangular, since the latter is essentially canonical. Type
\bprog
  M = H * qflll(H)
@eprog
Its column give an LLL-reduced basis for $\Lambda$ (\kbd{qflll(H)} itself
gives the base change matrix). The LLL algorithm outputs a nice basis for a
lattice given by an arbitrary basis, where nice means the basis vectors are
almost orthogonal and short, with precise guarantees on their relations to
the shortest vectors. Not really spectacular on this example, though.

Let us try something else, there should be an integer relation between
$\log 3$, $\log 5$ and $\log 15$. How to detect it?
\bprog
  u = [log(15), log(5), log(3)];
  m = matid(3); m[3,] = round(u * 10^25);
  v = qflll(m)[,1] \\@com first vector of the LLL-reduced basis
  u * v
@eprog\noindent
Pretty close. In fact, \kbd{lindep} automates this kind of search for integer
relations, using more sophisticated algorithms like PSLQ.

Let us come back to $\Lambda$ above, and our LLL basis in \kbd{M}. Type
\bprog
  G = M~*M  \\@com Gram matrix
  m = qfminim(G, norml2(M[,1]), 100);
@eprog\noindent
This enumerates the vectors in $\Lambda$ which are shorter than the first LLL
basis vector, at most 100 of them. There are $\kbd{m[1]} = 6$ such vectors,
and \kbd{m[3]} gives half of them (\kbd{-m[3]} would complete the lot): they
are the first 3 basis vectors! So these are optimally short, at least with
respect to the Euclidean length. Let us try
\bprog
  m = qfminim(G, norml2(M[,4]), 100, 2);
@eprog\noindent
(The flag $2$ instructs \kbd{qfminim} to use a different enumeration
strategy, which is much faster when we expect more short vectors than we want
to store. Without the flag, this example requires several hours. This is an
exponential time algorithm, after all!) This time, we find a slew of short
vectors; \kbd{matrank(m[3])} says the 100 we have are all included in a
2-dimensional space. Let us try
\bprog
  m = qfminim(G, norml2(M[,4]) - 1, 100000, 2);
@eprog\noindent
This time we find 50840 vectors of the requested length, spanning a
$4$-dimensional space, which is actually generated by \kbd{M[,1]},
\kbd{M[,2]} \kbd{M[,3]} and \kbd{M[,5]}.

\section{Using Transcendental Functions}

All the elementary transcendental functions and several higher transcendental
functions are provided: $\Gamma$ function, incomplete $\Gamma$ function, error
function, exponential integral, Bessel functions ($H^1$, $H^2$, $I$, $J$,
$K$, $N$), confluent hypergeometric functions, Riemann $\zeta$ function,
polylogarithms, Weber functions, theta functions. More will be written if the
need arises.

In this type of functions, the default precision plays an essential role.
In almost all cases transcendental functions work in the following way.
If the argument is exact, the result is computed using the current
default precision. If the argument is not exact, the precision of the
argument is used for the computation. A note of warning however: even in this
case the \emph{printed} value is the current real format, usually the
same as the default precision. In the present chapter we assume that your
machine works with 32-bit long integers. If it is not the case, we leave it
to you as a good exercise to make the necessary modifications.

Let's assume that we have 28 decimals of default precision (this is what we
get automatically at the start of a \kbd{gp} session on 32-bit machines). Type
\kbd{e = exp(1)}. We get the number $e=2.718\dots$ to 28 decimals. Let us check
how many correct decimals we really have. Change the precision to a
substantially higher value, for example by typing \kbd{\b{p} 50}. Then type
\kbd{e}, then \kbd{exp(1)} once again. This last value is the correct value
of the mathematical constant $e$ to 50 decimals, while the variable \kbd{e}
shows the value that was computed to 28 decimals. Clearly they coincide to
exactly 29 significant digits.

So 28 digits are printed, but how many significant digits are actually
contained in the variable \kbd{e}? Type \kbd{\#e} which indicates we have
exactly $3$ mantissa words. Since $3\ln(2^{32}) / \ln(10)\approx28.9$ we see
that we have 28 or 29 significant digits (on 32-bit machines).

\smallskip
Come back to 28 decimals (\kbd{\b{p} 28}). If we type \kbd{exp(1.)}
you can check that we also obtain 28 decimals. However, type
\kbd{f = exp(1 + 1E-30)}. Although the default precision is still 28,
you can check using the method above that we have in fact 59 significant
digits! The reason is that \kbd{1 + 1E-30} is computed according to the PARI
philosophy, i.e.~to the best possible precision. Since \kbd{1E-30} has 29
significant digits and 1 has ``infinite'' precision, the number \kbd{1 +
1E-30} will have $59=29+30$ significant digits, hence \kbd{f} also.

Now type \kbd{cos(1E-15)}. The result is printed as $1.0000\dots$, but
is of course not exactly equal to 1. Using \kbd{\#\%}, we see that the
result has 7 mantissa words, giving us the possibility of having 67
correct significant digits. In fact (look in precision 100), only 60 are
correct. PARI gives you as much as it can, and since 6 mantissa words
would have given you only 57 digits, it uses 7. But why does it give so
precise a result? Well, it is the same reason as before. When $x$ is close
to 1, $\cos(x)$ is close to $1-x^2/2$, hence the precision is going to be
approximately the same as when computing this quantity, here $1-0.5*10^{-30}$
where $0.5*10^{-30}$ is considered with 28 significant digit accuracy. Hence the
result will have approximately $28+30=58$ significant digits.

Unfortunately, this philosophy cannot go too far. For example, when you
type \kbd{cos(0)}, \kbd{gp} should give you exactly 1. Since it is reasonable for
a program to assume that a transcendental function never gives you an exact
result, \kbd{gp} gives you $1.000\dots$ with as many mantissa word as the current
precision.
\medskip
Let's see some more transcendental functions at work. Type
\kbd{gamma(10)}. No problem (type \kbd{9!} to check). Type \kbd{gamma(100)}.
The number is now written in exponential format because the default
accuracy is too small to give the correct result (type \kbd{99!} to get the
complete answer).
To get the complete value, there are two solutions. The first and most natural
one is to increase the precision. Since \kbd{gamma(100)} has 156 decimal
digits, type \kbd{\b{p} 170} to be on the safe side, then \kbd{gamma(100)}
once again. After some work, the printed result is this time perfectly
correct.

However, this is probably not the best way to proceed. Come back first to the
default precision (type \kbd{\b{p} 28}). As the gamma function increases
very rapidly, one usually uses its logarithm. Type \kbd{lngamma(100)}. This
time the result has a reasonable size, and is exactly equal to \kbd{log(99!)}.

Try \kbd{gamma(1/2 + 10*I)}. No problem, we have the complex gamma function.
Now type
\bprog
  t = 1000;
  z = gamma(1 + I*t) * t^(-1/2) * exp(Pi/2*t) / sqrt(2*Pi)
  norm(z)
@eprog\noindent The latter is very close to 1, in accordance with the complex
Stirling formula.
\smallskip

Let's play now with the Riemann zeta function. First turn on the timer (type
\kbd{\#}). Type \kbd{zeta(2)}, then \kbd{Pi\pow 2/6}. This seems correct. Type
\kbd{zeta(3)}. All this takes essentially no time at all. However, type
\kbd{zeta(3.)}. Although the result is the same, you will notice that the
time is substantially larger; if your machine is too fast to see the
difference, increase the precision to \kbd{\b{p}1000}. This is because PARI
uses special formulas to compute \kbd{zeta(n)} when \kbd{n} is an integer.

Type \kbd{zeta(1 + I)}. This also works. Now for fun, let us compute in a
very naive way the first complex zero of \kbd{zeta}. We know that it is
of the form $1/2 + i*t$ with $t$ between 14 and 15. Thus, we can use the
following series of instructions. But instead of typing them directly, write
them into a file, say \kbd{zeta.gp}, then type \kbd{\b{r} zeta.gp} under
\kbd{gp} to read it in:
\bprog
  {
    t1 = 1/2 + 14*I;
    t2 = 1/2 + 15*I; eps = 1E-50;
    z1 = zeta(t1);
    until (norm(z2) < eps,
      z2 = zeta(t2);
      if (norm(z2) < norm(z1),
        t3 = t1; t1 = t2; t2 = t3; z1 = z2
      );
      t2 = (t1+t2) / 2.;
      print(t1 ": " z1)
    )
  }
@eprog\noindent
Don't forget the braces: they tell \kbd{gp} that a sequence of instructions
is going to span many lines. We thus obtain the first zero to 25 significant
digits.

By the way, you don't need to type in the suffix~\kbd{.gp} in the \b{r}
command: it is supplied by \kbd{gp} if you forget it. The suffix is not
mandatory either, but it is convenient to have all GP scripts labeled in the
same distinctive way. Also, some text editors, e.g. Emacs or Vim, will
recognize GP scripts as such by their suffix and load special colourful
modes. \medskip
%
As mentioned at the beginning of this tutorial, some transcendental functions
can also be applied to $p$-adic numbers. This is as good a time as any to
familiarize yourself with them. Type
\bprog
  a = exp(7 + O(7^10))
  log(a)
@eprog\noindent All seems in order.
\bprog
  b = log(5 + O(7^10))
  exp(b)
@eprog\noindent Is something wrong? We don't recover the number we started
with? This is normal: type
\bprog
  exp(b) * teichmuller(5 + O(7^10))
@eprog\noindent
and we indeed recover our initial number. The Teichm\"uller
character \kbd{teichmuller(x)} on $\Z_p^*$ is the unique \hbox{$(p-1)$-st}
root of unity which is congruent to \kbd{x} modulo $p$, assuming that \kbd{x}
is a $p$-adic unit.\smallskip
%
Let us come back to real numbers for the moment. Type \kbd{agm(1,sqrt(2))}.
This gives the arithmetic-geometric mean of 1 and $\sqrt2$, and is the basic
method for computing complete elliptic integrals. In fact, type

\kbd{Pi/2 / intnum(t=0,Pi/2, 1 / sqrt(1 + sin(t)\pow 2))},

\noindent and the result is the same. The elementary transformation
\kbd{x = sin(t)} gives the mathematical equality
$$\int_0^1 \dfrac{dx}{\sqrt{1-x^4}} = \dfrac{\pi}{2\text{AGM}(1,\sqrt2)}
\enspace,$$
which was one of Gauss's remarkable discoveries in his youth.

Now type \kbd{2 * agm(1,I) / (1+I)}. As you see, the complex AGM also works,
although one must be careful with its definition. The result found is
almost identical to the previous one. Do you see why?

Finally, type \kbd{agm(1, 1 + 7 + O(7\pow 10))}. So we also have $p$-adic
AGM. Note however that since the square root of a $p$-adic number is not
in general an element of the same $p$-adic field,
only certain $p$-adic AGMs can be computed. In addition,
when $p=2$, the congruence restriction is that \kbd{agm(a,b)} can be computed
only when \kbd{a/b} is congruent to 1 modulo $16$, and not 8 as could be
expected.\smallskip
%
Now type \kbd{?3}. This gives you the list of all transcendental functions.
Instead of continuing with more examples, we suggest that you experiment
yourself with the list of functions. In each case, try integer, real, complex
and $p$-adic arguments. You will notice that some have not been implemented
(or do not have a reasonable definition).

\section{Using Numerical Tools}

 Although not written to be a numerical analysis package, PARI can
nonetheless perform some numerical computations. Since linear algebra and
polynomial computations are treated somewhere else, this section focuses on
solving equations and various methods of summation.

You of course know the formula $\pi = 4(1-\dfrac13+\dfrac15-\dfrac17+\cdots)$
which is deduced from the power series expansion of \kbd{atan(x)}. You also
know that $\pi$ cannot be computed from this formula, since the convergence
is so slow. Right? Wrong! Type
\bprog
  \p 100
  4 * sumalt(k=0, (-1)^k/(2*k + 1))
@eprog\noindent In a split second, we get $\pi$ to 100 significant digits
(type \kbd{Pi} to check).

Similarly, try
\bprog
  sumpos(k=1, k^-2)
@eprog\noindent Although once again the convergence is slow, the summation is
rather fast; compare with the exact result \kbd{Pi\pow 2/6}. This is less
impressive because a bit slower than for alternating sums, but still useful.

Even better, \kbd{sumalt} can be used to sum divergent series! Type
\bprog
  zet(s) = sumalt(k=1, (-1)^(k-1) / k^s) / (1 - 2^(1-s))
@eprog\noindent
Then for positive values of \kbd{s} different from 1, \kbd{zet(s)} is equal
to \kbd{zeta(s)} and the series converges, albeit slowly; \kbd{sumalt}
doesn't care however. For negative \kbd{s}, the series diverges, but
\kbd{zet(s)} still gives the correct result! (Namely, the value of a suitable
analytic continuation.) Try \kbd{zet(-1)}, \kbd{zet(-2)}, \kbd{zet(-1.5)},
and compare with the corresponding values of \kbd{zeta}. You should not push
the game too far: \kbd{zet(-100)}, for example, gives a completely wrong
answer.

Try \kbd{zet(I)}, and compare with \kbd{zeta(I)}. Even (some) complex values
work, although the sum is not alternating any more! Similarly, try
\bprog
  sumalt(n=1, (-1)^n / (n+I))
@eprog
\medskip
More traditional functions are the numerical integration functions.
Come back to \kbd{\b{p} 28} since these routines are very slow when working
with too many significant digits. Try \kbd{intnum(t=1,2, 1/t)}
and presto! you get 26 decimals of $\log(2)$. Look at Chapter 3 to see the
available integration functions.

With PARI, however, you can go further since complex types are allowed.
For example, assume that we want to know the location of the zeros of the
function $h(z)=e^z-z$. We use Cauchy's theorem, which tells us that the
number of zeros in a disk of radius $r$ centered around the origin is
equal to
$$\dfrac{1}{2i\pi}\int_{C_r}\dfrac{h'(z)}{h(z)}\,dz\enspace,$$
where $C_r$ is the circle of radius $r$ centered at the origin.
The function we want to integrate is
\bprog
  fun(z) = my(u = exp(z)); (u-1) / (u-z)
@eprog\noindent
(Here \kbd{u} is a local variable to the function \kbd{f}: whenever
a function is called, \kbd{gp} fills its argument list with the actual arguments
given, and initializes the other declared parameters and local variables to
0. It will then restore their former values upon exit. If we had not declared
\kbd{u} in the function prototype, it would be considered as a global
variable, whose value would be permanently changed. It is not mandatory to
declare in this way all parameters, but beware of side effects!)

Type now:
\bprog
  zero(r) = r/(2*Pi) * intnum(t=0, 2*Pi, real( fun(r*exp(I*t)) * exp(I*t) ))
@eprog
The function \kbd{zero(r)} will count the number of zeros of \kbd{fun}
whose modulus is less than \kbd{r}: we simply made the change of variable
$z = r*\exp(i*t)$, and took the real part to avoid integrating the
imaginary part. Actually, there is a built-in function \kbd{intcirc}
to integrate over a circle, yielding the much simpler:
\bprog
  zero(r) = intcirc(z=0, r, fun(z))
@eprog
(This is a little faster than the previous implementation, and no less
accurate.)

We may type \kbd{\b{p} 9} since we know that the result is a small integer
(but the computations should be instantaneous even at \b{p} 100 or so),
then \kbd{zero(1)}, \kbd{zero(1.5)}. The result tells us that there are
no zeros inside the unit disk, but that there are two (necessarily
complex conjugate) in the annulus $1<|z|<1.5$. For the sake of
completeness, let us compute them. Let $z = x+iy$ be such a zero, with
$x$ and $y$ real. Then the equation $e^z-z=0$ implies, after elementary
transformations, that $e^{2x}=x^2+y^2$ and that $e^x\cos(y)=x$. Hence
$y=\pm\sqrt{e^{2x}-x^2}$ and hence $e^x\cos(\sqrt{e^{2x}-x^2})=x$.
Therefore, type
\bprog
  fun(x) = my(u = exp(x)); u * cos(sqrt(u^2 - x^2)) - x
@eprog\noindent
Then \kbd{fun(0)} is positive while \kbd{fun(1)} is negative. Come
back to precision 28 and type
\bprog
  x0 = solve(x=0,1, fun(x))
  z = x0 + I*sqrt(exp(2*x0) - x0^2)
@eprog\noindent
which is the required zero. As a check, type \kbd{exp(z) - z}.

Of course you can integrate over contours which are more complicated than
circles, but you must perform yourself the changes of variable as we have
done above to reduce the integral to a number of integrals on line segments.
\smallskip
%
The example above also shows the use of the \kbd{solve} function. To use
\kbd{solve} on functions of a complex variable, it is necessary to reduce the
problem to a real one. For example, to find the first complex zero of the
Riemann zeta function as above, we could try typing

\kbd{solve(t=14,15, real( zeta(1/2 + I*t) ))},

\noindent but this does not work because the real part is positive for
$\kbd{t}=14$ and $15$. As it happens, the imaginary part works. Type

\kbd{solve(t=14,15, imag( zeta(1/2 + I*t) ))},

\noindent and this now works. We could also narrow the search interval and
type for instance

\kbd{solve(t=14,14.2, real( zeta(1/2 + I*t) ))}

\noindent which would also work.

\section{Polynomials}

First a word of warning: it is essential to understand the difference between
exact and inexact objects. Try
\bprog
  gcd(x - Pi, x^2 - 6*zeta(2))
@eprog\noindent
We return a trivial GCD because the notion of GCD for non-exact polynomials
doesn't make much sense. A better quantitative approach is to use
\bprog
  polresultant(x - Pi, x^2 - 6*zeta(2))
@eprog\noindent
A result close to zero shows that the GCD is non-trivial for small
deformations of the inputs. Without telling us what it is, of course. This
being said, we will mostly use polynomials (and power series) with exact
coefficients in our examples.\smallskip

The simplest way to input a polynomial, is to simply write it down,
or use an explicit formula for the coefficients and the function \kbd{sum}:
\bprog
  T = 1 + x^2 + 27*x^10;
  T = sum(i = 1, 100, (i+1) * x^i);
@eprog\noindent but it is in much more efficient to create a vector of
coefficients then convert it to a polynomial using \kbd{Pol} or \kbd{Polrev}
(\kbd{Pol([1,2])} is $x+2$, \kbd{Polrev([1,2]) is $2x + 1$}) :
\bprog
  T = Polrev( vector(100, i, i) );

  for (i=1, 10^4, Polrev( vector(100, i, i) ) )   \\@com time: 124ms
  for (i=1, 10^4, sum(i = 1, 100, (i+1) * x^i) ) \\@com time: 3,985ms
@eprog\noindent The reason for the discrepancy is that the explicit summation
(of densely encoded polynomials) is quadratic in the degree, whereas creating
a vector of coefficients then converting it to a polynomial type is linear.

We also have a few built-in classical polynomial families. Consider the
$15$-th cyclotomic polynomial,
\bprog
  pol = polcyclo(15)
@eprog\noindent
which is of degree $\varphi(15)=8$. Now, type
\bprog
  r = polroots(pol)
@eprog\noindent We obtain the 8 complex roots of \kbd{pol}, given to 28
significant digits. To see them better, type \b{b}: they are given as pairs
of complex conjugate roots, in a random order. The only ordering done by the
function \kbd{polroots} concerns the real roots, which are given first, and
in increasing order.

The roots of \kbd{pol} are by definition the primitive $15$-th roots of unity.
To check this, simply type \kbd{rc = r\pow 15}. Why, we get an error message!
Fair enough, vectors cannot be multiplied, even less raised to a power.
However, type
\bprog
  rc = r^15.
@eprog\noindent without forgetting the `\kbd{.}' at the end. Now it works,
because powering to a non-integer exponent is a transcendental function and
hence is applied termwise. Note that the fact that $15.$ is a real number
which is representable exactly as an integer has nothing to do with the
problem.

We see that the components of the result are very close to 1. It is however
tedious to look at all these real and imaginary parts. It would be impossible
if we had many more. Let's do it automatically. Type
\bprog
  rr = round(rc)
  sqrt( norml2(rc - rr) )
@eprog\noindent We see that \kbd{rr} is indeed all 1's, and that the
$L^2$-norm of \kbd{rc - rr} is around $3.10^{-28}$, reasonable enough when we
work with 28 significant digits! Note that the function \kbd{norml2},
contrary to what its name implies, does not give the $L^2$ norm but its
square, hence we must take the square root. Well, this is not absolutely
necessary in the present case! In fact, \kbd{round} itself already provides
a built-in rough approximation of the error:
\bprog
  rr = round(rc, &e)
@eprog\noindent Now \kbd{e} contains the number of error bits when rounding
\kbd{rc} to \kbd{rr}; in other words the sup norm of $\kbd{rc} - \kbd{rr}$
is bounded by $2^{-\kbd{e}}$.
%
\smallskip
Now type
\bprog
  pol = x^5 + x^4 + 2*x^3 - 2*x^2 - 4*x - 3
  factor(pol)
  factor( poldisc(pol) )
  fun(p) = factorpadic(pol,p,10);
@eprog\noindent The polynomial \kbd{pol} factors over $\Q$ (or $\Z$) as a
product of two factors, and the primes dividing its discriminant are
$11$, $23$ and $37$. We also created a function \kbd{fun(p)} which factors
\kbd{pol} over $\Q_p$ to $p$-adic precision 10. Type
\bprog
  fun(5)
  fun(11)
  fun(23)
  fun(37)
@eprog\noindent to see different splittings.

Similarly, type
\bprog
  lf(p) = lift( factormod(pol,p) );
  lf(2)
  lf(11)
  lf(23)
  lf(37)
@eprog\noindent which show the different factorizations, this time over
$\F_p$. In fact, even better: type successively
\bprog
  t = 't;   \\@com we want \kbd{t} to be a free variable
  T = ffinit(3,3, t)
  pol2 = subst(T, t, x)
  fq = factorff(pol2, 3, T)
  centerlift( lift(fq) )
@eprog\noindent
\kbd{T}, which is actually \kbd{t\pow 3 + t\pow 2 + t + 2} (with intmod
coefficients), is defined above to be an irreducible polynomial of degree $3$
over $\F_3$. This code snippet factors the polynomial \kbd{pol2} over the
finite field $\F_3[t]/(T)$. This is of course a form of the field $\F_{27}$.
We know that Gal$(\F_{27}/\F_3)$ is cyclic of order 3 generated by the
Frobenius homomorphism $u\mapsto u^3$, and the roots that we have found give
the action of the powers of the Frobenius on \kbd{t}. (If you do not know
what I am talking about, please try some examples, it's not so hard to figure
out.) We took pain above to factor a polynomial in the variable $x$ over a
finite field defined by a polynomial in $t$, even though they were apparently
one and the same. There is a crucial rule in all routines involving relative
extensions: the variable associated to the base field is required to have
lower priority than the variables of polynomials whose coefficients are taken
in that base field. Have a look at the section on \emph{Variable priorities}
in the user's manual (see ``The GP programming language'').

Similarly, type
\bprog
  pol3 = x^4 - 4*x^2 + 16
  fn = lift( factornf(pol3, t^2 + 1) )
@eprog\noindent and we get the factorization of the polynomial \kbd{pol3}
over the number field defined by $t^2+1$, i.e.~over $\Q(i)$. Without the
\kbd{lift}, the result would involve number field elements as \typ{POLMOD}s
of the form \kbd{Mod(1+t, t\pow2+1)}, which are more explicit but much less
readable.
\smallskip

To summarize, in addition to being able to factor integers, you can
factor polynomials over $\C$ and $\R$ using \kbd{polroots},
over $\F_p$ using \kbd{factormod}, over $\F_{p^k}$ using \kbd{factorff},
over $\Q_p$ using \kbd{factorpadic}, over $\Q$ using \kbd{factor}, and over
number fields using \kbd{factornf} or \kbd{nffactor}. Note however
that \kbd{factor} itself will try to guess intelligently over which ring you
want to factor: try
\bprog
  pol = x^2 + 1;
  factor(pol)
  factor(pol *1.)
  factor(pol * (1 + 0*I))
  factor(pol * (1 + 0.*I))
  factor(pol * Mod(1,2))
  factor(pol * Mod(1, Mod(1,3)*(t^2+1)))
@eprog

In the present version \vers{}, it is not possible to factor over other
rings than the ones mentioned above. For example \kbd{gp} cannot directly
factor multivariate polynomials, although it is not too hard to write a
simple-minded Kronecker's substitution to reduce to the univariate case.
(Exercise!) Other functions related to factoring are \kbd{padicappr},
\kbd{polrootsmod}, \kbd{polrootspadic}, \kbd{polsturm}. Play with them a
little.

Finally, type
\bprog
  polsym(pol3, 20)
@eprog\noindent where \kbd{pol3} was defined above. This gives the sum of the
$k$-th powers of the roots of \kbd{pol3} up to $k=20$, of course computed
using Newton's formula and not using \kbd{polroots}. You notice that every
odd sum is zero (expected, since the polynomial is even), but also that
the signs follow a regular pattern and that the (non-zero) absolute values
are powers of 2. This is true: prove it, and more precisely find an explicit
formula for the $k$-th symmetric power not involving (non-rational) algebraic
numbers.

\section{Power series}

Now let's play with power series as we have done at the beginning.  Type
\bprog
  N = 39;
  8*x + prod(n=1,N, if(n%4, 1 - x^n, 1), 1 + O(x^(N+1)))^8
@eprog\noindent
Apparently, only even powers of \kbd{x} appear. This is surprising, but can
be proved using the theory of modular forms. Note that we initialize
the product to \kbd{1 + O(x\pow (N+1))}, otherwise the whole computation would
be done with polynomials; this would first have been slightly slower and also
totally useless since the coefficients of \kbd{x\pow (N+1)} and above are
irrelevant anyhow if we stop the product at $\kbd{n} = \kbd{N}$.

While we are on the subject of modular forms (which, together with Taylor
series expansions are another great source of power series), type
\bprog
  \ps 122     \\@com shortcut for \kbd{default(seriesprecision, 122)}
  d = x * eta(x)^24
@eprog\noindent This gives the first 122 terms of the Fourier series
expansion of the modular discriminant function $\Delta$ of Ramanujan. Its
coefficients give by definition the Ramanujan $\tau$ function, which has a
number of marvelous properties (look at any book on modular forms for
explanations). We would like to see its properties modulo 2. Type \kbd{d\%2}.
Hmm, apparently PARI doesn't like to reduce coefficients of power series, or
polynomials for that matter, directly. Can we do it without writing a little
program? No problem. Type instead
\bprog
  lift(Mod(1,2) * d)
  centerlift(Mod(1,3) * d)
@eprog\noindent and now this works like a charm. The pattern in the first
result is clear; the pattern is less clear in the second result, but
nonetheless there is one. Of course, it now remains to prove it (see Antwerp
III or your resident modular forms guru).

\section{Working with Elliptic Curves}

Now we are getting to more complicated objects. Just as with number fields
which we will meet later on, the first thing to do is to initialize them.
That's because a lot of data will be needed repeatedly, and it's much more
convenient to have it ready once and for all. Here, this is done with the
function \kbd{ellinit}.

So type
\bprog
  e0 = ellinit([6,-3,9,-16,-14])
@eprog
This computes a number of things
about the elliptic curve defined by the affine equation
%
$$ y^2+6xy+9y = x^3-3x^2-16x-14\enspace. $$
%
It is not that clear what all these funny numbers mean, except that we
recognize the first few of them as the coefficients we just input. To
retrieve meaningful information from such complicated objects (and number
fields will be much worse), one uses so-called \emph{member
functions}. Type \kbd{?.} to get a complete list. Whenever \kbd{ell} appears
in the right hand side, we can apply the corresponding function to an object
output by \kbd{ellinit}. (I'm sure you know how the other \kbd{init} functions
will be called now, don't you? Oh, by the way, neither \kbd{clgpinit} nor
\kbd{pridinit} exist.)

  Let's try it. The discriminant \kbd{e0.disc} is equal to 37, hence
the conductor of the curve is 37. Of course in general it is not so
trivial. In fact, although the equation of the curve is clearly
minimal (since the discriminant is $12$th-power-free), it is not in
standard reduced form, so type
\bprog
  e = ellminimalmodel(e0)
@eprog\noindent which
gives the \kbd{ell} structure associated with the standard model,
exactly as if we had used \kbd{ellinit} on a reduced equation. For
some related data, type
\bprog
  gr = ellglobalred(e0)
@eprog\noindent The first
component \kbd{gr[1]} tells us that the conductor is 37 as we already
knew.  The second component is a 4-component vector which allows us to
get the minimal equation: in fact \kbd{e} is \kbd{ellchangecurve(e0, gr[2])}.
You can for the moment ignore the third component \kbd{gr[3]}. (For the
impatient reader, this is the product of the local Tamagawa numbers, $c_p$.)
Type
\bprog
  q0 = [-2,2]
  ellisoncurve(e0, q0)
  q = ellchangepoint(q0,gr[2])
  ellisoncurve(e, q)
@eprog\noindent
The point \kbd{q0} is on the curve, as checked by \kbd{ellisoncurve}, and we
transferred it onto the minimal model \kbd{e}, using \kbd{ellchangepoint} and
the change of variable computed above. Note that \kbd{ellchangepoint()} is
unusual among the elliptic curve functions in that it does not take an
\kbd{ell} structure as its first argument: in \kbd{gp}, points do not
``know'' which curve they are on, but to move a point from one model to
another we only need to know the coordinates of the point and the
transformation data here stored in \kbd{gr[2]}.  Also, the point at infinity
is represented as \kbd{[0]} on all elliptic curves; this is the identity for
the group law.

Here, \kbd{q=[0,0]} obviously lies on \kbd{e}, which has equation
$y^2+y = x^3-x$.  Let us now play a little with points on \kbd{e}.
The group law on an elliptic curve is implemented with the functions
\kbd{elladd} for addition, \kbd{ellsub} for subtraction and
\kbd{ellpow} for multiplication by an integer (``powering'').  For
example, the negative of \kbd{q} is \kbd{ellsub(e,[0],q)}, and the
double is obtained either as \kbd{ellpow(e,q,2)} or as
\kbd{elladd(e,q,q)}.

Now \kbd{q} may be a torsion point. Type \kbd{ellheight(e, q)}, which
computes the canonical Neron-Tate height of \kbd{q}. Note that
\kbd{ellheight} assumes that \kbd{e} is \emph{minimal}! (Which it is.)
This is non-zero, hence \kbd{q} is not torsion. To see this even better,
type
\bprog
  for(k = 1, 20, print(ellpow(e, q, k)))
@eprog\noindent and we see the characteristic parabolic explosion of the size
of the points. (And another proof that \kbd{q} is not torsion, assuming
Mazur's bound on the size of the rational torsion.) We could can also type
\kbd{ellorder(e, q)} which returns 0, telling us yet again that \kbd{q} is
non-torsion. As a consistency check, type
\bprog
  ellheight(e, ellpow(e, q,20)) / ellheight(e, q)
@eprog\noindent
We indeed find $400=20^2$ as it should be.

Notice how (almost) all those \kbd{ell}--prefixed functions take our
elliptic curve as a first argument? This will be true with number
fields as well: whatever object was initialized by an $ob$--\kbd{init}
function will have to be used as a first argument of all the
$ob$--prefixed functions. Conversely, you won't be able to use any
such high-level function before you correctly initialize the relevant
object. \smallskip

Ok, let's try another curve. Type
\bprog
  E = ellinit([0,-1,1,0,0])
  q = [0,0]; ellheight(E, q)
@eprog\noindent This corresponds to the equation $y^2+y = x^3-x^2$ and an
obvious rational point on it. Again from the discriminant we see that the
conductor is equal to 11, and if you type \kbd{ellminimalmodel(E)} you will
see that the equation  for \kbd{E} is minimal. This time the height is very
close to zero, hence \kbd{q} must be a torsion point. Indeed, typing
\bprog
  for(k=1, 5, print(ellpow(E, q,k)))
  ellorder(E, q)   \\@com simpler
@eprog\noindent
we see in two different ways that \kbd{q} is a point of order 5. Moreover,
typing
\bprog
  elltors(E)
@eprog\noindent shows that \kbd{q} generates all the torsion of \kbd{E},
which is cyclic of order~$5$. \smallskip

Let's try still another curve, $y^2+y = x^3-7x+6$:
\bprog
  e = ellinit([0,0,1,-7,6])
  ellglobalred(e)
@eprog\noindent As before, this is a minimal equation; now the conductor is
5077. There are some trivial integral points on this curve, but let's try to
be more systematic. Typing
\bprog
  elltors(e)
@eprog\noindent
shows that the torsion subgroup is trivial, so we don't have to worry about
torsion points. Next, the member \kbd{e.roots} gives us the 3 roots of the
minimal equation over $\C$, i.e.~$Y^2=X^3-7X+25/4$ (set $(X,Y)=(x,y+1/2)$) so
if $(x,y)$ is a real point on the curve, $x$ must be at least equal to the
smallest root, i.e.~$x\ge-3$. Finally, if $(x,y)$ is on the curve, its
opposite is clearly $(x,-y-1)$. Type
\bprog
  {
    v = List();
    for (x = -3, 1000,
      s = ellordinate(e,x);
      if (#s, listput(v, [x,s[1]]))  \\@com or \kbd{if (\#s != 0,}
    ); v
  }
@eprog\noindent
(If cardinality of \kbd{s} is non-zero, insert the first point with given
\kbd{x}.) We thus get a large number (18) of integral points. Together with
their opposites and the point at infinity, this makes a total of 37 integral
points, which is large for a curve having such a small conductor. So we
suspect (if we do not know already, since this curve is quite famous!) that
the rank of this curve must be large. Let's try and put some order into this.
Note that we work only with the integral points, but in general rational
points should also be considered. Type
\bprog
  v = Vec(v);  \\@com convert to a vector
  hv = ellheight(e, v)
@eprog\noindent (We convert the list to a vector because \typ{LIST}s are not
accepted as input by most mathematical functions. They must be converted to a
standard vector type first.) This gives the vector of canonical heights. Let
us order the points according to their height. For this, type
\bprog
  perm = vecsort(vector(#v,i,i), (a,b) -> hv[a]>hv[b]);
  v = vecextract(v, perm)
@eprog\noindent
Note how we input a non-obvious comparison function to \kbd{vecsort}, using
an inline anonymous function: we sort the integers between $1$ and the number
of points we have declaring that $a > b$ if and only if the height of the
$a$-th point is larger than the height of the $b$-th. The result is a
permutation which we apply to \kbd{v} using \kbd{vecextract}.

It seems reasonable to take the numbers with smallest height as
possible generators of the Mordell-Weil group. Let's try the first
four: type
\bprog
  m = ellheightmatrix(e, vecextract(v,[1,2,3,4])); matdet(m)
@eprog\noindent
Since the curve has no torsion, the determinant being close to zero
implies that the first four points are dependent. To find the
dependency, it is enough to find the kernel of the matrix \kbd{m}. So
type \kbd{matker(m)}: we indeed get a non-trivial kernel, and the
coefficients are close to integers. Typing \kbd{elladd(e, v[1],v[3])} does
indeed show that it is equal to \kbd{v[4]}.

Taking any other four points, we would in fact always find a
dependency.  Let's find all dependencies. Type
\bprog
  vp = [v[1],v[2],v[3]]~;
  m = ellheightmatrix(e,vp);
  matdet(m)
@eprog\noindent This is now clearly non-zero so the first 3 points are
linearly independent, showing that the rank of the curve is at least equal to
3. (It is in fact equal to 3, and \kbd{e} is the curve of smallest conductor
having rank 3.) We would like to see whether the other points are dependent.
For this, we use the function \kbd{ellbil}. Indeed, if \kbd{Q} is some point
which is dependent on \kbd{v[1],v[2]} and \kbd{v[3]}, then \kbd{matsolve(m,
ellbil(e, vp,Q))} will by definition give the coefficients of the dependence
relation. If these coefficients are close to integers, then there is a
dependency, otherwise not.  This is much safer than using the \kbd{matker}
function. Thus, type
\bprog
  w = vector(#v, k, matsolve(m, ellbil(e, vp,v[k])))
  wr = round(w, &err)
@eprog\noindent
We ``see'' that the coefficients are all very close to integers, and we
quantified it with the last instruction: \kbd{wr} is the vector expressing
all the components of \kbd{v} on its first 3, and \kbd{err} gives an upper
bound on the maximum distance to an integer. We are thus led to strongly
believe that the curve has rank exactly 3, with generators \kbd{v[1]},
\kbd{v[2]} and \kbd{v[3]}, and this can be proved to be the case.

Two remarks: (1) Using the height pairing to find dependence relations
as we have done only in fact finds relations modulo torsion; but in
this example, the curve has trivial torsion, as we checked
earlier. (2) In the above calculation we were lucky that all the
\kbd{v[j]} were $\Z$-linear combinations of \kbd{v[1]}, \kbd{v[2]} and
\kbd{v[3]} and not just $\Q$-linear combinations;  in general the
result of \kbd{matsolve(m, ellbil(e, vp,Q))} might have given a vector
of rationals: if $k\ge2$ is minimal such that $kQ$ is in the subgroup
generated by \kbd{v[1]}, \kbd{v[2]} and \kbd{v[3]}, then the entries of
\kbd{matsolve(m, ellbil(e, vp,Q))} will be rationals with common
denominator~$k$.

\smallskip

Let us explore a few more elliptic curve related functions. Keep our
curve \kbd{e} of rank 3, and type
\bprog
  v1 = [1,0]; z1 = ellpointtoz(e, v1)
  v2 = [2,0]; z2 = ellpointtoz(e, v2)
@eprog\noindent
We thus get the complex parametrization of the curve. To add the points
\kbd{v1} and \kbd{v2}, we should of course type \kbd{elladd(e, v1,v2)},
but we can also type \kbd{ellztopoint(e, z1 + z2)} which has the disadvantage
of giving complex numbers, but illustrates how the group law on \kbd{e} is
obtained from the addition law on $\C$.

Type
\bprog
  f = x * Ser(ellan(e, 30))
@eprog\noindent This gives a power series which is the Fourier expansion of a
modular form of weight 2 for $\Gamma_0(5077)$. (This has been proved
directly, before Wiles's general result.) In fact, to find the modular
parametrization of the curve, type
\bprog
  modul = elltaniyama(e)
  u = modul[1];
  v = modul[2];
@eprog\noindent
We can check in various ways that this indeed parametrizes the curve:
\bprog
  (v^2 + v) - (u^3 - 7*u + 6)
@eprog\noindent
is close to $0$ for instance, or simply that \kbd{ellisoncurve(e,modul)}
returns~$1$. Now type
\bprog
  x * u' / (2*v + 1)
@eprog\noindent and we see that this is equal to the modular form \kbd{f}
found above; the quote \kbd{'} tells \kbd{gp} to take the derivative of the
expression with respect to its main variable. The functions \kbd{u} and
\kbd{v}, considered on the upper half plane with $x=e^{2i\pi\tau}$, are in
fact modular \emph{functions} for $\Gamma_0(5077)$. \smallskip

The function \kbd{ellan(e, 30)} gives the first~$30$ coefficients
$a_n$ of the $L$-series of \kbd{e}.  One can also ask for a single
coefficient: the millionth is \kbd{ellak(e, 10\pow 6)}.  Note however
that calling \kbd{ellan(e,100000)} is much faster than
the equivalent \kbd{vector(100000,k,ellak(e,k))}.  For a
prime~\kbd{p},
\kbd{ellap(e,p)} is equivalent to \kbd{ellak(e,p)};  this is the
integer $a_p$ such that the number of points on \kbd{e} over $\F_p$ is
$1+p-a_p$. (With the standard PARI distribution, \kbd{ellap} is the only way
to obtain the order of an elliptic curve over $\F_p$ in \kbd{gp}. The
external package \kbd{ellsea} provides much more efficient routines.)

Finally, let us come back to the curve \kbd{E} defined above by
\kbd{E = ellinit([0,-1,1,0,0])}, which had an obvious rational $5$-torsion
point. The sign of the functional equation, given by \kbd{ellrootno(E)}, is
$+1$. Assuming the rank parity conjecture, it follows that the Mordell-Weil
group of $E$ has even rank. The value of the L-function of \kbd{E} at 1
\bprog
  ls = elllseries(E, 1)
@eprog\noindent is definitely non-zero, so \kbd{E} has rank $0$. According to
the Birch and Swinnerton-Dyer conjecture, which is proved for this curve,
\kbd{ls} is given by the following formula (in this case):
%
\def\sha{\hbox{III}}
$$L(E,1)=\dfrac{\Omega\cdot c\cdot|\sha|}{|E_{\text{tors}}|^2}\enspace,$$
%
where $\Omega$ is the real period of $E$, $c$ is the global Tamagawa number,
product of the local $c_p$ for primes $p$ dividing the conductor, $|\sha|$ is
the order of the Tate-Shafarevich group, and $E_{\text{tors}}$ is the
torsion group of $E$.

Now we know many of these quantities: $\Omega$ is equal to \kbd{E.omega[1]}
(if there had been 3 real roots instead of 1 in \kbd{E.roots}, $\Omega$ would
be equal to \kbd{2 * E.omega[1]}). The Tamagawa number $c$ is given as the
last component of \kbd{ellglobalred(E)}, and here is equal to 1. We already
know that the torsion subgroup of $E$ contains a point of order 5, and typing
\kbd{elltors(E)} shows that it is of order exactly 5. So type
\bprog
  ls * 25/E.omega[1]
@eprog\noindent This shows that $\sha$ must be the trivial group.

For more detailed information on the local reduction of an elliptic curve at
a specific prime~\kbd{p}, use the function \kbd{elllocalred(E,p)}; the second
component gives the Kodaira symbol in an encoded form.  See the manual or
online help for details.

\section{Working in Quadratic Number Fields}

The simplest of all number fields outside $\Q$ are quadratic fields and are
the subject of the present section. We shall deal in the next one with
general number fields (including $\Q$ and quadratic fields!), and one should
be aware that all we will see now has a more powerful, in general easier to
use, equivalent in the general context. But possibly much slower.

Such fields are characterized by their discriminant. Even better, any
non-square integer $D$ congruent to 0 or 1 modulo 4 is the discriminant of a
specific order in a quadratic field. We can check whether this order is
maximal with \kbd{isfundamental(D)}. Elements of a quadratic field are of the
form $a+b\omega$, where $\omega$ is chosen as $\sqrt{D}/2$ if $D$ is even and
$(1+\sqrt{D})/2$ if $D$ is odd, and are represented in PARI by quadratic
numbers. To initialize working in a quadratic order, one starts by the
command \kbd{w = quadgen($D$)}.

This sets \kbd{w} equal to $\omega$ as above, and is printed \kbd{w}. Note
however that if several different quadratic orders are used, a printed \kbd{w}
may have several different meanings. For example if you type
\bprog
  w1 = quadgen(-23)
  w2 = quadgen(-15)
@eprog\noindent
both will be printed as \kbd{w}, but of course they are not equal. Hence
beware when dealing with several quadratic orders at once. \smallskip
%
In addition to elements of a quadratic order, we also want to be able to
handle ideals of such orders. In the quadratic case, it is equivalent to
handling binary quadratic forms. For negative discriminants, quadratic forms
are triples $(a,b,c)$ representing the form $ax^2+bxy+cy^2$. Such a form will
be printed as, and can be created by, \kbd{Qfb($a$,$b$,$c$)}.

Such forms can be multiplied, divided, powered as many PARI objects using
the usual operations, and they can also be reduced using the function
\kbd{qfbred} (it is not the purpose of this tutorial to explain what all
these things mean). In addition, Shanks's NUCOMP algorithm has been
implemented (functions \kbd{qfbnucomp} and \kbd{qfbnupow}), and this is
usually a little faster.

Finally, you have at your disposal the functions \kbd{qfbclassno} which
(\emph{usually}) gives the class number, the function \kbd{qfbhclassno}
which gives the Hurwitz class number, and the much more sophisticated
\kbd{quadclassunit} function which gives the class number and class group
structure.

Let us see examples of all this at work.

Type \kbd{qfbclassno(-10007)}. \kbd{gp} tells us that the result is 77. However,
you may have noticed in the explanation above that the result is only
\emph{usually} correct. This is because the implementers of the algorithm
have been lazy and have not put the complete Shanks algorithm into PARI,
causing it to fail in certain very rare cases. In practice, it is almost
always correct, and the much more powerful \kbd{quadclassunit} program, which
\emph{is} complete (at least for fundamental discriminants) can give
confirmation; but now, under the Generalized Riemann Hypothesis!

So we may be a little suspicious of this class number. Let us check it.
First, we need to find a quadratic form of discriminant $-10007$. Since this
discriminant is congruent to 1 modulo 8, we know that there is an ideal of
norm equal to 2, i.e.~a binary quadratic form $(a,b,c)$ with $a=2$. To
compute it we type \kbd{f = qfbprimeform(-10007, 2)}. OK, now we have a form.
If the class number is correct, the very least is that this form raised to
the power 77 should equal the identity. Type \kbd{f\pow 77}. We get a form
starting with 1, i.e.~the identity. Raising \kbd{f} to the powers 11 and 7
does not give the identity, thus we now know that the order of \kbd{f} is
exactly 77, hence the class number is a multiple of 77. But how can we be
sure that it is exactly 77 and not a proper multiple? Well, type
\bprog
  sqrt(10007)/Pi * prodeuler(p=2,500, 1./(1 - kronecker(-10007,p)/p))
@eprog\noindent
This is nothing else than an approximation to the Dirichlet class number
formula. The function \kbd{kronecker} is the Kronecker symbol, in this case
simply the Legendre symbol. Note also that we have written \kbd{1./(1 - \dots)}
with a dot after the first 1. Otherwise, PARI may want to compute the whole
thing as a rational number, which would be terribly long and useless. In fact
PARI does no such thing in this particular case (\kbd{prodeuler} is always
computed as a real number), but you never know. Better safe than sorry!

We find 77.77, pretty close to 77, so things seem in order. Explicit bounds
on the prime limit to be used in the Euler product can be given which make
the above reasoning rigorous.

Let us try the same thing with $D=-3299$. \kbd{qfbclassno} and the Euler
product convince us that the class number must be 27. However, we get stuck
when we try to prove this in the simple-minded way above. Indeed, we type
\kbd{f = qfbprimeform(-3299, 3)} (2 is not the norm of a prime ideal but
3 is), and we see that \kbd{f} raised to the power 9 is equal to the identity.
This is the case for any other quadratic form we choose. So we suspect that the
class group is not cyclic. Indeed, if we list all 9 distinct powers of \kbd{f},
we see that \kbd{qfbprimeform(-3299, 5)} is not on the list, although its cube
is as it must. This implies that the class group is probably equal to a
product of a cyclic group of order 9 by a cyclic group of order 3. The Euler
product plus explicit bounds prove this.

Another way to check it is to use the \kbd{quadclassunit} function by typing
for example
\bprog
  quadclassunit(-3299)
@eprog\noindent
Note that this function cheats a little and could still give a wrong answer,
even assuming GRH: the forms given could generate a strict subgroup of the
class group. If we want to use proven bounds under GRH, we have to type
\bprog
  quadclassunit(-3299,,[1,6])
@eprog\noindent
The double comma \kbd{,,} is not a typo, it means we omit an optional second
argument. As we want to use the optional \emph{third} argument, we have to
indicate to \kbd{gp} we skipped this one.

Now, if we believe in GRH, the class group is as we thought (see Chapter 3
for a complete description of this function).

  Note that using the even more general function \kbd{bnfinit} (which handles
general number fields and gives more complicated results), we could
\emph{certify} this result, i.e~remove the GRH assumption. Let's do it, type
\bprog
  bnf = bnfinit(x^2 + 3299); bnfcertify(bnf)
@eprog

  A non-zero result (here 1) means that everything is ok. Good, but what did
we certify after all? Let's have a look at this \kbd{bnf} (just type it!).
Enlightening, isn't it? Recall that the \kbd{init} functions (we've already
seen \kbd{ellinit}) store all kind of technical information which you
certainly don't care about, but which will be put to good use by higher level
functions. That's why \kbd{bnfcertify} could not be used on the output of
\kbd{quadclassunit}: it needs much more data.

  To extract sensible information from such complicated objects, you must use
one of the many \emph{member functions} (remember: \kbd{?.} to get a complete
list). In this case \kbd{bnf.clgp} which extracts the class group structure.
This is much better. Type \kbd{\%.no} to check that this leading 27 is indeed
what we think it is and not some stupid technical parameter. Note that
\kbd{bnf.clgp.no} would work just as well, or even \kbd{bnf.no}!

As a last check, we can request a relative equation for the Hilbert class
field of $\Q(\sqrt{-3299})$: type \kbd{quadhilbert(-3299)}. It is indeed of
degree 27 so everything fits together.

\medskip
%
Working in real quadratic fields instead of complex ones, i.e.~with $D>0$, is
not very different.

The same \kbd{quadgen} function is used to create elements. Ideals are again
represented by binary quadratic forms $(a,b,c)$, this time indefinite. However,
the Archimedean valuations of the number field start to come into play, hence
in fact quadratic forms with positive discriminant will be represented as a
quadruplet $(a,b,c,d)$ where the quadratic form itself is $ax^2+bxy+cy^2$
with $a$, $b$ and $c$ integral, and $d\in \R$ is a ``distance''
component, as defined by Shanks and Lenstra.

To create such forms, one uses the same function as for definite ones, but
you can add a fourth (optional) argument to initialize the distance:
\kbd{q = Qfb($a$, $b$, $c$, $d$)}.
If the discriminant of \kbd{poldisc(q)} is negative, $d$ is silently
discarded. If you omit it, this component is set to \kbd{0.} (i.e.~a real
zero to the current precision).

Again these forms can be multiplied, divided, powered, and they can be
reduced using \kbd{qfbred}. This function is in fact a succession of
elementary reduction steps corresponding essentially to a continued fraction
expansion, and a single one of these steps can be achieved by adding an
(optional) flag to the arguments of \kbd{qfbred}. Since handling the
fourth component $d$ usually involves computing expensive logarithms, the
same flag may be used to ignore the fourth component. Finally, it is
sometimes useful to operate on forms of positive discriminant without
performing any reduction (this is useless in the negative case), the
functions \kbd{qfbcompraw} and \kbd{qfbpowraw} do exactly that.

Again, the function \kbd{qfbprimeform} gives a prime form, but the form which
is given corresponds to an ideal of prime norm which is usually not reduced.
If desired, it can be reduced using \kbd{qfbred}.

Finally, you still have at your disposal the function \kbd{qfbclassno} which
gives the class number (this time \emph{guaranteed} correct),
\kbd{quadregulator} which gives the regulator, and the much more sophisticated
\kbd{quadclassunit} giving the class group's structure and its generators,
as well as the regulator. The \kbd{qfbclassno} and \kbd{quadregulator}
functions use an algorithm which is $O(\sqrt D)$, hence become very slow for
discriminants of more than 10 digits. \kbd{quadclassunit} can be used on a
much larger range.

Let us see examples of all this at work and learn some little known number
theory at the same time. First of all, type
\bprog
  d = 3 * 3299; qfbclassno(d)
@eprog\noindent We see that the class number is 3. (We know
in advance that it must be divisible by 3 from the \kbd{d = -3299} case above
and Scholz's mirror theorem.) Let us create a form by typing
\bprog
  f = qfbred(qfbprimeform(d,2), 2)
@eprog\noindent
(the last 2 tells \kbd{qfbred} to ignore the Archimedean component). This
gives us a prime ideal of norm equal to 2. Is this ideal principal? Well, one
way to check this, which is not the most efficient but will suffice for now,
is to look at the complete cycle of reduced forms equivalent to \kbd{f}. Type
\bprog
 g = f; for(i=1,20, g = qfbred(g, 3); print(g))
@eprog\noindent
(this time the 3 means to do a single reduction step, still not using
Shanks's distance). We see that we come back to the form \kbd{f} without
having the principal form (starting with $\pm1$) in the cycle, so the ideal
corresponding to \kbd{f} is not principal.

Since the class number is equal to 3, we know however that \kbd{f\pow 3} will
be a principal ideal $\alpha\Z_K$. How do we find $\alpha$? For this, type
\bprog
  f3 = qfbpowraw(f, 3)
@eprog
This computes the cube of \kbd{f}, without
reducing it. Hence it corresponds to an ideal of norm equal to $8=2^3$, so we
already know that the norm of $\alpha$ is equal to $\pm8$. We need more
information, and this will be given by the fourth component of the form.
Reduce your form until you reach the unit form (you will have to type
\kbd{qfbred(\%,~1)} exactly 6 times), then extract the Archimedean component,
say $c$. By definition of this distance, we know that
$${\alpha\over{\sigma(\alpha)}}=\pm e^{2c},$$
where $\sigma$ denotes real conjugation in our quadratic field. This can be
automated:
\bprog
  q = f3;
  while(abs(component(q,1)) != 1, print(q); q = qfbred(q, 1))
  c = component(q,4);
@eprog\noindent
Thus, if we type
\bprog
  a = sqrt(8 * exp(2*c))
  sa = 8 / a
@eprog\noindent
we know that up to sign, \kbd{a} and \kbd{sa} are
numerical approximations of $\alpha$ and $\sigma(\alpha)$. Of course,
$\alpha$ can always be chosen to be positive, and a quick numerical check
shows that the difference of \kbd{a} and \kbd{sa} is close to an integer, and
not the sum, so that in fact the norm of $\alpha$ is equal to $-8$ and the
numerical approximation to $\sigma(\alpha)$ is \kbd{$-$sa}. Thus we type
\bprog
  p = x^2 - round(a-sa)*x - 8
@eprog\noindent
and this is the characteristic polynomial of $\alpha$. We can check that the
discriminant of this polynomial is a square multiple of \kbd{d}, so $\alpha$
is indeed in our field. More precisely, solving for $\alpha$ and using the
numerical approximation that we have to resolve the sign ambiguity in the
square root, we get explicitly $\alpha=(15221+153\sqrt d)/2$. Note that this
can also be done automatically using the functions \kbd{polred} and
\kbd{modreverse}, as we will see later in the general number field case, or
by solving a system of 2 linear equations in 2 variables. (Exercise: now that
we have $\alpha$, check that it is indeed a generator of the ideal
corresponding to the form \kbd{f3}.)

\medskip Let us now play a little with cycles. Type
\bprog
  D = 10^7 + 1
  quadclassunit(D)
@eprog\noindent
We get as a result a 5-component vector, which tells us that (under GRH) the
class number is equal to 1, and the regulator is approximately
equal to $2641.5$. You may certify this with
\bprog
  bnf = bnfinit(x^2 - D, 1);  \\ @com insist on finding fundamental unit
  bnfcertify(bnf);
@eprog\noindent
although it's a little inefficient. Indeed \kbd{bnfcertify} needs the
fundamental unit which is so large that \kbd{bnfinit} will have a hard time
computing it: it needs about $R/\log(10)\approx 1147$ digits of precision!
(So that it would have given up had we not inserted the flag $1$.)
See \kbd{bnf.fu}. On the other hand, you can try \kbd{quadunit(D)}.
Impressive, isn't it? (You can check that its logarithm is indeed equal to
the regulator.)

Now just as an example, let's assume that we want the regulator to 500
decimals, say. (Without cheating and computing the fundamental unit exactly
first!) I claim that simply from the crude approximation above, this can be
computed with no effort.

This time, we want to start with the unit form. Type:
\bprog
  u = qfbred(qfbprimeform(D, 1), 2)
@eprog\noindent
We use the function \kbd{qfbred} with no distance since we want the initial
distance to be equal to~0. Now type  \kbd{f = qfbred(u, 1)}. This is the
first form encountered along the principal cycle. For the moment, keep the
precision low, for example the initial default precision. The distance from
the identity of \kbd{f} is around 4.253. Very crudely, since we want a
distance of $2641.5$, this should be encountered approximately at
$2641.5/4.253=621$ times the distance of \kbd{f}. Hence, as a first try, we
type \kbd{f\pow 621}. Oops, we overshot, since the distance is now $3173.02$.
Now we can refine our initial estimate and believe that we should be close to
the correct distance if we raise \kbd{f} to the power $621*2641.5/3173$ which
is close to $517$. Now if we compute \kbd{f\pow 517} we hit the principal
form right on the dot. Note that this is not a lucky accident: we will always
land extremely close to the correct target using this method, and usually at
most one reduction correction step is necessary. Of course, only the distance
component can tell us where we are along the cycle.

Up to now, we have only worked to low precision. The goal was to obtain this
unknown integer $517$. Note that this number has absolutely no mathematical
significance: indeed the notion of reduction of a form with positive
discriminant is not well defined since there are usually many reduced forms
equivalent to a given form. However, when PARI makes its computations, the
specific order and reductions that it performs are dictated entirely by the
coefficients of the quadratic form itself, and not by the distance component,
hence the precision used has no effect.

Hence we now start again by setting the precision to (for example) 500,
we retype the definition of \kbd{u} (why is this necessary?), and then
\kbd{qfbred(u, 1)\pow 517}. Of course we know in advance that we land on the
unit form, and the fourth component gives us the regulator to 500 decimal
places with no effort at all.

In a similar way, we could obtain the so-called \emph{compact representation}
of the fundamental unit itself, or $p$-adic regulators. I leave this as
exercises for the interested reader.

You can try the \kbd{quadhilbert} function on that field but, since the class
number is $1$, the result won't be that exciting. If you try it on our
preceding example ($3*3299$) it should take about 2 seconds.
\medskip

Time for a coffee break?

\section{Working in General Number Fields}
\subsec{Elements}

The situation here is of course more difficult. First of all, remembering
what we did with elliptic curves, we need to initialize data linked to our
base number field, with something more serious than \kbd{quadgen}. For
example assume that we want to work in the number field $K$ defined by one of
the roots of the equation $x^4+24x^2+585x+1791=0$. This is done by typing
\bprog
  T = x^4 + 24*x^2 + 585*x + 1791
  nf = nfinit(T)
@eprog\noindent
We get an \kbd{nf} structure but, thanks to member functions, we do not need
to know anything about it. If you type \kbd{nf.pol}, you will get the
polynomial \kbd{T} which you just input. \kbd{nf.sign} yields the signature
$(r_1,r_2)$ of the field, \kbd{nf.disc} the field discriminant, \kbd{nf.zk}
an integral basis, etc\dots.

The integral basis is expressed in terms of a generic root \kbd{x} of \kbd{T}
and we notice it's very far from being a power integral basis, which is a
little strange for such a small field. Hum, let's check that: \kbd{poldisc(T)}?
Ooops, small wonder we had such denominators, the index of the power order
$\Z[x]/(T)$ in the maximal order $\Z_K$ is, well,
\kbd{nf.index}. Note that this is also given by
\bprog
  sqrtint(poldisc(nf.pol) / nf.disc)
@eprog\noindent
Anyway, that's $3087$, we don't want to work with such a badly skewed
polynomial! So, we type
\bprog
  P = polred(T)
@eprog\noindent
We see from the third component that the
polynomial $x^4-x^3-21x^2+17x+133$ defines the same field with much smaller
coefficients, so type
\bprog
  A = P[3]
@eprog\noindent
The \kbd{polred} function usually gives a simpler polynomial, and also
sometimes some information on the existence of subfields. For example in this
case, the second component of \kbd{polred} tells us that the field defined by
$x^2-x+1=0$, i.e.~the field generated by the cube roots of unity, is a subfield
of~$K$. Note this is incidental information and that the list of subfields
found in this way is usually far from complete. To get the complete list, one
uses \kbd{nfsubfields} (we shall do that later on).

  Type \kbd{poldisc(A)}, this is much better, but maybe not optimal yet
(the index is still $7$). Type \kbd{polredabs(A)} (the \kbd{abs} stands for
absolute). Since it seems that we won't get anything better, we'll stick with
\kbd{A} (note however that \kbd{polredabs} finds a smallest generating
polynomial with respect to a bizarre norm which ensures that the index will
be small, but not necessarily minimal). In fact, had you typed
\kbd{nfinit(T, 3)}, \kbd{nfinit} would first have tried to find a good
polynomial defining the same field (i.e.~one with small index) before
proceeding.

  It's not too late, let's redefine our number field:
\bprog
  NF = nfinit(nf, 3)
@eprog\noindent
The output is a two-component vector. The first component is the new
\kbd{nf}: type
\bprog
  nf = NF[1];
@eprog\noindent
If you type \kbd{nf.pol}, you notice that \kbd{gp} indeed replaced your bad
polynomial \kbd{T} by a much better one, which happens to be \kbd{A}. (Small
wonder, \kbd{nfinit} internally called \kbd{polredabs}.) The second component
enables you to switch conveniently to our new polynomial.

Namely, call $\theta$ a root of our initial polynomial \kbd{T}, and $\alpha$
a root of the one that \kbd{polred} has found, namely \kbd{A}. These are
algebraic numbers, and as already mentioned are represented as polmods. For
example, in our special case $\theta$ and $\alpha$ are equal to the polmods
\bprog
  THETA = Mod(x, x^4 + 24*x^2 + 585*x + 1791)
  ALPHA = Mod(x, x^4 - x^3 - 21*x^2 + 17*x + 133)
@eprog\noindent respectively. Here we are considering only the algebraic
aspect, and hence ignore completely \emph{which} root $\theta$ or $\alpha$ is
chosen.

Now you may have a number of elements of your number field which are
expressed as polmods with respect to your old polynomial, i.e.~as polynomials
in $\theta$. Since we are now going to work with $\alpha$ instead, it is
necessary to convert these numbers to a representation using $\alpha$. This
is what the second component of \kbd{NF} is for: type \kbd{C = NF[2]}, you get
\bprog
  Mod(-10/7*x^3 + 43/7*x^2 + 73/7*x - 64, x^4 - x^3 - 21*x^2 + 17*x + 133)
@eprog\noindent
meaning that $\theta =
-\dfrac{10}{7}\alpha^3+\dfrac{43}{7}\alpha^2+\dfrac{73}{7}\alpha-64$, and hence the conversion from a
polynomial in $\theta$ to one in $\alpha$ is easy, using \kbd{subst}. (We
could get this polynomial from \kbd{polred} as well, try \kbd{polred(T, 2)}.)
If we want the reverse, i.e.~to go back from a representation in $\alpha$ to
a representation in $\theta$, we use the function \kbd{modreverse} on this
polmod \kbd{C}. Try it. The result has a big denominator (1029) essentially
because our initial polynomial \kbd{T} was so bad. By the way, to get that
1029,
you should type \kbd{denominator(content(C))}. Trying \kbd{denominator} by
itself would not work since the denominator of a polynomial is defined to be
1 (and its numerator is itself). The reason for this is that we think of a
polynomial as a special case of a rational function.\smallskip

From now on, we forget about \kbd{T}, and use only the polynomial
\kbd{A} defining $\alpha$, and the components of the vector \kbd{nf} which
gives information on our number field $K$. Type
\bprog
  u = Mod(x^3 - 5*x^2 - 8*x + 56, A) / 7
@eprog\noindent
This is an element in $K$. There are three equivalent
representations for number field elements: polmod, polynomial, and column
vector giving a decomposition in the integral basis \kbd{nf.zk} (\emph{not} on
the power basis $(1,x,x^2,\dots)$). All three are equally valid when the
number field is understood (is given as first argument to the function).
You will be able to use any one of them as long as the function you call
requires an \kbd{nf} argument as well. However, most PARI functions will
return elements as column vectors. It's an important feature of number
theoretic functions that, although they may have a preferred format for
input, they will accept a wealth of other different formats. We already saw
this for \kbd{nfinit} which accepts either a polynomial or an \kbd{nf}. It
will be true for ideals, congruence subgroups, etc.

  Let's stick with elements for the time being. How does one go from one
representation to the other? Between polynomials and polmods, it's easy:
\kbd{lift} and \kbd{Mod} will do the job. Next, from polmods/polynomials to
column vectors: type \kbd{v = nfalgtobasis(nf, u)}. So $\kbd{u} = \alpha^3-
\alpha^2 - \alpha + 8$, right? Wrong! The coordinates of \kbd{u} are given
with respect to the \emph{integral basis}, not the power basis
$(1,\alpha,\alpha^2,\alpha^3)$ (and they don't coincide, type \kbd{nf.zk} if
you forgot what the integral basis looked like). As a polynomial in $\alpha$,
we simply have $\kbd{u} = {1\over7}(\alpha^3 - 5\alpha^2-8\alpha+56)$, which
is trivially deduced from the original polmod representation!

Of course \kbd{v = nfalgtobasis(nf, lift(u))} would work equally well. Indeed
we don't need the polmod information since \kbd{nf} already provides the
defining polynomial. To go back to polmod representation, use
\kbd{nfbasistoalg(nf, v)}. Notice that \kbd{u} is an algebraic integer since
\kbd{v} has integer coordinates (try \kbd{denominator(v) == 1}, which is of
course overkill here, but not so in a program).

Let's try this out. We may for instance compute \kbd{u\pow 3}. Try it. Or, we
can type \kbd{1/u}. Better yet, if we want to know the norm from $K$ to $\Q$
of \kbd{u}, we type \kbd{norm(u)} (what else?); \kbd{trace(u)} works as well.
Notice that none of this would work on polynomials or column vectors since
you don't have the opportunity to supply \kbd{nf}! But we could use
\kbd{nfeltpow(nf,u,3)}, \kbd{nfeltdiv(nf,1,u)} (or \kbd{nfeltpow(nf,u,-1)})
which would work whatever representation was chosen. There is no
\kbd{nfeltnorm} function (\kbd{nfelttrace} does not exist either), but we can
easily write one:
\bprog
  nfeltnorm(nf,u) = norm(nfbasistoalg(nf, u))
@eprog\noindent
Notice that this is certainly not foolproof (try it with complex or quadratic
arguments!), but we could refine it if the need arose. You can also
consider $(u)$ as a principal ideal, and just type
\bprog
  idealnorm(nf,u)
@eprog\noindent
Of course, in this way, we lose the \emph{sign} information. We will talk
about ideals later on.

  If we want all the symmetric functions of \kbd{u} and not only the norm, we
type \kbd{charpoly(u)}. Note that this gives the characteristic polynomial of
\kbd{u}, and not in general the minimal polynomial.
\misctitle{Exercise} How does one get the minimal polynomial from this?
Find a simpler expression for \kbd{u}. \smallskip

  Now let's work on the field itself. The \kbd{nfinit} command already
gave us some information. The field is totally complex (its signature
\kbd{nf.sign} is $[0,2]$), its discriminant \kbd{nf.disc} is $18981$ and
\kbd{nf.zk} is an integral basis. The Galois group of its Galois closure
can be obtained by typing \kbd{polgalois(A)}. The answer (\kbd{[8,-1,1]}, or
\kbd{[8,-1,1,"D(4)"]} if the \kbd{galdata} package is installed) shows that
it is equal to $D_4$, the dihedral group with 8 elements, i.e.~the group of
symmetries of a square.

This implies that the field is ``partially Galois'', i.e.~that there exists
at least one non-trivial field isomorphism which fixes $K$, exactly one in
this case. Type \kbd{nfgaloisconj(nf)}. The result tells us that, apart from
the trivial automorphism, the map
$$\alpha \mapsto {1\over7}(-\alpha^3+5\alpha^2+\alpha-49)$$
is the only field automorphism.
\bprog
  nfgaloisconj(nf);
  s = Mod(%[2], A)
  charpoly(s)
@eprog\noindent
and we obtain \kbd{A} once again.

Let us check that \kbd{s} is of order 2: \kbd{subst(lift(s), x, s)}. It is. We
may express it as a matrix:
\bprog
  w = Vec( matid(4) ) \\@com canonical basis
  v = vector(#w, i, nfgaloisapply(nf, s, w[i]))
  M = Mat(v)
@eprog\noindent
The vector \kbd{v} contains the images of the integral basis elements (as
column vectors). The last statement concatenates them into a square matrix.
So, \kbd{M} gives the action of \kbd{s} on the integral basis. Let's check
\kbd{M\pow2}. That's the identity all right.

The fixed field of this automorphism is going to be the only non-trivial
subfield of $K$. I seem to recall that \kbd{polred} told us this was the
third cyclotomic field. Let's check this: type \kbd{nfsubfields(nf)}. Indeed,
there's a quadratic subfield, but it's given by \kbd{T = x\pow 2 + 22*x + 133
} and I don't recognize it. But \kbd{nfisisom(T, polcyclo(3))} indeed tells
us that the fields $\Q[x]/(T)$ and $\Q[x]/(x^2+x+1)$ are isomorphic.
(In fact, \kbd{polred(T)} would tell us the same, but does not correspond to
a foolproof test: \kbd{polred} could have returned some other polynomials.)

We may also check that \kbd{k = matker(M-1)} is two-dimensional, then \kbd{z
= nfbasistoalg(nf, k[,2])} generates the quadratic subfield. Notice that 1,
\kbd{z} and \kbd{u} are $\Q$-linearly dependent, and in fact $\Z$-linearly as
well. Exercise: how would you check these two assertions in general? (Answer:
\kbd{concat}, then respectively \kbd{matrank} or \kbd{matkerint} (or
\kbd{qflll})). \kbd{z = charpoly(z)}, \kbd{z = gcd(z,z')} and \kbd{polred(z)}
tell us that we found back the same subfield again (as we ought to!).

Final check: type \kbd{nfrootsof1(nf)}. Again we find that $K$ contains
a cube root of unity, since the torsion subgroup of its unit group
has order 6. The given generator happens to be equal to \kbd{u}.

\misctitle{Additional comment} (you're no longer supposed to skip this,
but do as you wish):

Before working with ideals, let us note one more thing. The main part of the
work of \kbd{polred} or \kbd{nfinit} is to compute an integral basis, i.e.~a
$\Z$-basis of the maximal order $\Z_K$ of $K$. For a large polynomial, this
implies factoring the discriminant of the polynomial, which is often not
feasible. The situation may be improved in many ways:

1) First, it is often the case that our number field is of quite a special
type. For example, one may know in advance some prime divisors of the
discriminant. Hence we can ``help'' PARI by giving it that information. More
precisely, we can use the function \kbd{addprimes} to inform PARI to keep on
eye for these prime numbers. Do it only for big primes (bigger than
\kbd{primelimit}, whose value you can get using \kbd{default})~--- it is
useless otherwise.

2) The second way in which the situation may be improved is that often we do
not need the complete information on the maximal order, but only require that
the order be $p$-maximal for a certain number of primes $p$ (but then, we
may not be able to use functions which require a genuine \kbd{nf}). The
function \kbd{nfbasis} specifically computes the integral basis and is not
much quicker than \kbd{nfinit} so is not very useful in its standard use. But
you can provide a factorization of the discriminant as an optional third
argument. And here we can cheat, and give on purpose an incomplete
factorization involving only the primes we want. For example coming back to
our initial polynomial $T$, the discriminant of the polynomial is
$3^7\cdot7^6\cdot19\cdot37$. If we only want a $7$-maximal order, we simply
type
\bprog
  nfbasis(T, ,[7,6; 1537461,1])
@eprog\noindent
and the factors of 1537461 will not be looked at! (In this example it would
be stupid to cheat, but if the discriminant has 2000 digits, this can be a
handy trick.)

3) A last way will work if the \emph{field} discriminant is smooth (never
mind the polynomial discriminant). Simply call \kbd{B = nfbasis(T, 1)}, which
will return a basis for an order which is possibly not maximal. In fact, it
will be maximal if the field discriminant is smooth with respect to the
precomputed prime table, including \kbd{addprimes} offsprings (assuming the
later are genuine primes). You may then input the resulting basis to
\kbd{nfinit}, as \kbd{nfinit([T, B])}.

\subsec{Ideals}

We now want to work with ideals and not only
with elements. An ideal can be represented in many different ways. First, an
element of the field (in any of the various guises seen above) will be
considered as a principal ideal. Then the standard representation is a
square matrix giving the Hermite Normal Form (HNF) of a $\Z$-basis of the
ideal expressed on the integral basis \kbd{nf.zk}. Standard means that most
ideal related functions will use this representation for their output.

Prime ideals can be represented in a special form as well (see
\kbd{idealprimedec}) and all ideal-related functions will accept them. On the
other hand, the function \kbd{idealtwoelt} can be used to find a two-element
$\Z_K$-basis of a given ideal (as $a\Z_K + b\Z_K$, where $a$ and $b$ belong
to $K$), but this is \emph{not} a valid representation for an ideal under
\kbd{gp}, and most functions will choke on it (or worse, take it for
something else and output a meaningless result). To be able to use such an
ideal, you will first have to convert it to HNF form.

Whereas it's very easy to go to HNF form (use \kbd{idealhnf(nf,id)} for valid
ideals, or \kbd{idealhnf(nf,a,b)} for a two-element representation as above),
it's a much more complicated problem to check whether an ideal is principal
and find a generator. In fact an \kbd{nf} does not contain enough data for
this particular task. We'll need a Buchmann Number Field, or \kbd{bnf}, for
that. In particular, we need the class group and fundamental units, at least
in some approximate form. More on this later (which will trivialize the end
of the present section).\smallskip

Let us keep our number field $K$ as above and its \kbd{nf} structure. Type
\bprog
  P = idealprimedec(nf,7)
@eprog\noindent
This gives the decomposition of the prime number 7 into prime ideals. We have
chosen 7 because it divides \kbd{nf.index} (in fact, is equal to it), hence
is the most difficult case to treat.

The result is a vector with 4 components, showing that 7 is totally split in
the field $K$ into prime ideals of norm 7 (you can check:
\kbd{idealnorm(nf,P[1])}). Let us take one of these ideals, say the first, so
type
\bprog
  pr = P[1]
@eprog
We obtain its inertia and residue degree as \kbd{pr.e}
and \kbd{pr.f}, and its two generators as \kbd{pr.gen}. One of them is
$\kbd{pr.p} = 7$, and the other is guaranteed to have valuation $1$ at
\kbd{pr}. What is the Hermite Normal Form of this ideal? No problem:
\bprog
  idealhnf(nf,pr)
@eprog\noindent and we have the desired HNF. Let's now perform ideal
operations. For example type
\bprog
  idealmul(nf, pr, idealmul(nf, pr,pr))
@eprog\noindent or more simply
\bprog
  pr3 = idealpow(nf, pr,3)
@eprog\noindent
to get the cube of the ideal \kbd{pr}. Since the norm of this ideal is equal
to $343=7^3$, to check that it is really the cube of \kbd{pr} and not of
other ideals above 7, we can type
\bprog
  for(i=1, #P, print( idealval(nf, pr3, P[i]) ))
@eprog\noindent
and we see that the valuation at \kbd{pr} is equal to 3, while the others are
equal to zero. We could see this as well from \kbd{idealfactor(nf, pr3)}.

Let us now work in the class group ``by hand'' (we shall see simpler ways
later). We shall work with \emph{extended ideals}: an extended ideal is a pair
\kbd{[A, t]}, where $A$ is an ordinary ideal as above, and $t$ a number field
element; this pair represents the ideal $(t) A$.
\bprog
  id3 = [pr3, 1]
  r0 = idealred(nf, id3)
@eprog\noindent
The input \kbd{id3} is an extended ideal: pr3 together with 1 (trivial
factorization). The new extended ideal \kbd{r0} is equal to the old one,
in the sense that the products $(t)A$ are the same. It contains a ``reduced''
ideal equivalent to \kbd{pr3} (modulo the principal ideals), and a generator
of the principal ideal that was factored out.

Now, just for fun type
\bprog
  r = r0; for(i=1,3, r = idealred(nf,r, [1,5]); print(r))
@eprog\noindent
The ideals in the third \kbd{r} and initial \kbd{r0} are equal, say
$(t) A = (t_0) A$: this means we
have found a unit $(t_0/t)$ in our field, and it is easy to extract this unit
given the extended component:
\bprog
  t0 = r0[2]; t = r[2];
  u = nfeltdiv(nf, t0, t)
  u = nfbasistoalg(nf, u)
@eprog\noindent The last line recovers the unit as an algebraic number.
Type
\bprog
  ch = charpoly(u)
@eprog\noindent
and we obtain the characteristic polynomial \kbd{ch} of $u$ again. (Whose
constant coefficient is $1$, hence $u$ is indeed a unit.)

There is of course no reason for $u$ to be a fundamental unit. Let us see if
it is a square. Type
\bprog
  F = factor(subst(ch, x, x^2))
@eprog\noindent
We see that \kbd{ch(x\pow2)} is a product of 2 polynomials of degree 4, hence
$u$ is a square. (Why?) We now want to find its square root. A simple method
is as follows:
\bprog
  NF = subst(nf,x,y);
  r = F[1,1] % (x^2 - nfbasistoalg(NF, u))
@eprog\noindent
to find the remainder of the characteristic polynomial of \kbd{u2} divided by
\kbd{x\pow 2 - $u$}. This is a polynomial of degree 1 in \kbd{x}, with polmod
coefficients, and we know that \kbd{u2}, being a root of both polynomials,
is the root of \kbd{r}, hence can be obtained by typing
\bprog
  u2 = -polcoeff(r,0) / polcoeff(r,1)
@eprog\noindent There is an important technicality in the above: why did we
need to substitute \kbd{NF} to \kbd{nf} ? The reason is that data related to
\kbd{nf} is given in terms of the variable \kbd{x}, seen modulo \kbd{nf.pol};
but we need \kbd{x} as a free variable for our polynomial divisions. Hence
the substitution of \kbd{x} by \kbd{y} in our \kbd{nf} data.

The most natural method is to try directly
\bprog
  nffactor(nf, y^2 - u)
@eprog\noindent
Except that this won't work for the same technical reason as above: the main
variable of the polynomial to be factored must have \emph{higher} priority
than the number field variable. This won't be possible here since \kbd{nf}
was defined using the variable \kbd{x} which has the highest possible
priority. So we need to substitute variables around:
\bprog
  nffactor(NF, x^2 - nfbasistoalg(NF, subst(lift(u),x,y)))
@eprog\noindent
(Of course, with better planning, we would just have defined \kbd{nf} in
terms of the \kbd{'y} variable, to avoid all these substitutions.)
\smallskip

A much simpler approach is to consider the above as instances of a
\emph{discrete logarithm} problem, where we want to express some elements an
abelian group (of finite type) in terms of explicitly given generators, and
transfer all computations from abstract groups like $\text{Cl}(K)$ and
$\Z_K^*$ to products of simpler groups like $\Z^n$ or $\Z/d\Z$. We shall do
exactly that in the next section.

Before that, let us mention another famous (but in fact, simpler)
\emph{discrete logarithm} problem, namely the one associated to the
invertible elements modulo an ideal: $(\Z_K / I)^*$. Just use \kbd{idealstar}
(this is an \kbd{init} function) and \kbd{ideallog}.

Many more functions on ideals are available. We mention here the complete
list, referring to Chapter 3 for detailed explanations:

\kbd{idealadd}, \kbd{idealaddtoone}, \kbd{idealappr}, \kbd{idealchinese},
\kbd{idealcoprime}, \kbd{idealdiv}, \kbd{idealfactor}, \kbd{idealhnf},
\kbd{idealintersect}, \kbd{idealinv}, \kbd{ideallist}, \kbd{ideallog},
\kbd{idealmin}, \kbd{idealmul}, \kbd{idealnorm}, \kbd{idealpow},
\kbd{idealprimedec}, \kbd{idealred},
\kbd{idealstar}, \kbd{idealtwoelt}, \kbd{idealval},
\kbd{nfisideal}.

We suggest you play with these to get a feel for the algebraic number theory
package. Remember that when a matrix (usually in HNF) is output, it is always
a $\Z$-basis of the result expressed on the \emph{integral basis} \kbd{nf.zk}
of the number field, which is usually \emph{not} a power basis.

\subsec{Class groups and units, \kbd{bnf}}

Apart from the above functions you have at your disposal the powerful
function \kbd{bnfinit}, which initializes a \kbd{bnf} structure, i.e.~a
number field with all its invariants (including class group and units), and
enough technical data to solve discrete logarithm problems in the class and
unit groups.

First type \kbd{setrand(1)}: this resets the random seed (to make sure we
and you get the exact same results). Now type
\bprog
  bnf = bnfinit(NF);
@eprog\noindent
where \kbd{NF} is the same number field as before. You do not want to see the
output clutter a number of screens so don't forget the semi-colon. (Well if
you insist, it is about three screenful in this case, but may require several
Megabytes for larger degrees.) Note that \kbd{NF} is now expressed in terms
of the variable \kbd{y}, to avoid later problems with variable priorities.

A word of warning: both the \kbd{bnf} and all results obtained from it are
\emph{conditional} on a Riemann Hypothesis at this point; the \kbd{bnf} must
be certified before the following statements become actual theorems.
\smallskip

Member functions are still available for \kbd{bnf} structures. So, let's try
them: \kbd{bnf.pol} gives \kbd{A}, \kbd{bnf.sign}, \kbd{bnf.disc},
\kbd{bnf.zk}, ok nothing really exciting. In fact, an \kbd{nf} is included
in the \kbd{bnf} structure: \kbd{bnf.nf} should be identical to
\kbd{NF}. Thus, all functions which took an \kbd{nf} as first argument, will
equally accept a \kbd{bnf} (and a \kbd{bnr} as well which contains even more
data).

Anyway, new members are available now: \kbd{bnf.no} tells us the class number
is 4, \kbd{bnf.cyc} that it is cyclic (of order 4 but that we already knew),
\kbd{bnf.gen} that it is generated by the ideal \kbd{g = bnf.gen[1]}. If you
\kbd{idealfactor(bnf, g)}, you recognize \kbd{P[2]}. (You may also play in
the other direction with \kbd{idealhnf}.) The regulator \kbd{bnf.reg} is
equal to $3.794\dots$. \kbd{bnf.tu} tells us that the roots of unity in $K$
are exactly the sixth roots of 1 and gives a primitive root $\zeta =
{1\over7}(\alpha^3 - 5\alpha^2 - 8\alpha + 56)$, which we have seen already.
Finally \kbd{bnf.fu} gives us a fundamental unit $\epsilon =
{1\over7}(\alpha^3 - 5\alpha^2 - \alpha + 28)$, which must be linked to the
units \kbd{u} and \kbd{u2} found above since the unit rank is~1. To find
these relations, type
\bprog
  bnfisunit(bnf, u)
  bnfisunit(bnf, u2)
@eprog\noindent
Lo and behold, \kbd{u = $\zeta^2\epsilon^2$} and \kbd{u2 =
$\zeta^{4}\epsilon^1$}.

\misctitle{Note} Since the fundamental unit obtained depends on the random
seed, you could have obtained another unit than $\epsilon$, had you not reset
the random seed before the computation. This was the purpose of the initial
\kbd{setrand} instruction, which was otherwise unnecessary.\medskip

We are now ready to perform operations in the class group. First and
foremost, let us certify the result: type \kbd{bnfcertify(bnf)}. The
output is \kbd{1} if all went well; in fact no other output is possible,
whether the input is correct or not, but you can get an error message (or in
exceedingly rare cases an infinite loop) if it is incorrect.

It means that we now know the class group and fundamental units
unconditionally (no more GRH then!). In this case, the certification process
takes a very short time, and you might wonder why it is not built in as a
final check in the \kbd{bnfinit} function. The answer is that as the
regulator gets bigger this process gets increasingly difficult, and becomes
soon impractical, while \kbd{bnfinit} still happily spits out results. So it
makes sense to dissociate the two: you can always check afterwards, if the
result is interesting enough. Looking at the tentative regulator, you know in
advance whether the certification can possibly succeed: if \kbd{bnf.reg} is
large, don't waste your time.


Now that we feel safe about the \kbd{bnf} output, let's do some real work.
For example, let us take again our prime ideal \kbd{pr} above 7. Since we
know that the class group is of order 4, we deduce that \kbd{pr} raised to
the fourth power must be principal. Type
\bprog
  pr4 = idealpow(nf, pr, 4)
  v = bnfisprincipal(bnf, pr4)
@eprog\noindent
The first component gives the factorization of the ideal in the class group.
Here, \kbd{[0]} means that it is up to equivalence equal to the 0-th power of
the generator \kbd{g} given in \kbd{bnf.gen}, in other words that it is a
principal ideal. The second component gives us the algebraic number $\alpha$
such that $\kbd{pr4}=\alpha\Z_K$, $\alpha$ being as usual expressed on the
integral basis. Type \kbd{alpha = v[2]}. Let us check that the result is
correct: first, type \kbd{idealnorm(bnf, alpha)}. (Note that we can use a
\kbd{bnf} with all the \kbd{nf} functions; but not the other way round, of
course.) It is indeed equal to $7^4 = 2401$, which is the norm of \kbd{pr4}.
This is only a first check. The complete check is obtained by computing the
HNF of the principal ideal generated by \kbd{alpha}. To do this, type
\kbd{idealhnf(bnf, alpha) == pr4}.

Since the equality is true, \kbd{alpha} is correct (not that there was any
doubt!). But \kbd{bnfisprincipal} also gives us information for non-principal
ideals. For example, type
\bprog
  v = bnfisprincipal(bnf, pr)
@eprog\noindent
The component \kbd{v[1]} is now equal to \kbd{[3]}, and tells us that \kbd{pr}
is ideal-equivalent to the cube of the generator \kbd{g}. Of course we
already knew this since the product of \kbd{P[3]} and \kbd{P[4]} was
principal (generated by \kbd{al}), as well as the product of all the
\kbd{P[$i$]} (generated by 7), and we noticed that \kbd{P[2]} was equal
to \kbd{g}, which has order 4. The second component \kbd{v[2]} gives us
$\alpha$ on the integral basis such that $\kbd{pr}=\alpha \kbd{g}^3$. Note
that if you \emph{don't} want this $\alpha$, which may be large and whose
computation may take some time, you can just add the flag $1$ (see the online
help) to the arguments of \kbd{bnfisprincipal}, so that it only returns the
position of \kbd{pr} in the class group. \smallskip

\subsec{Class field theory, \kbd{bnr}}

We now survey quickly some class field theoretic routines. We must first
initialize a Buchmann Number Ray, or \kbd{bnr}, structure, associated to a
\kbd{bnf} base field and a modulus. Let's keep $K$, and try a finite modulus
${\goth f} = 7\Z_K$. (See the manual for how to include infinite places in
the modulus.) Since $K$ will now become a base field over which we want to
build relative extensions, the associated \kbd{bnf} needs to have variables
of lower priority than the polynomials defining the extensions. Fortunately,
we already took care that, but it would have been easy to deal with the
problem now (as easy as \kbd{bnf = subst(bnf, x, y)}). Then type
\bprog
  bnr = bnrinit(bnf, 7, 1);
  bnr.cyc
@eprog\noindent
tells us the ray class group modulo ${\goth f}$ is isomorphic to
$\Z/24\Z \times \Z/6\Z \times \Z/2\Z $. The associated generators are
\kbd{bnr.gen}.  Just as a \kbd{bnf} contained an \kbd{nf}, a \kbd{bnr}
contains a \kbd{bnf} (hence an \kbd{nf}), namely \kbd{bnr.bnf}. Here
\kbd{bnr.clgp} refers to the ray class group, while \kbd{bnr.bnf.clgp} refers
to the class group.
\bprog
  rnfkummer(bnr,, 2)
  rnfkummer(bnr,, 3)
@eprog\noindent
outputs defining polynomials for the $2$ abelian extensions of $K$ of degree
$2$ (resp.~$3$), whose conductor is exactly equal to ${\goth f}$ (the modulus
used to define \kbd{bnr}). (In the current implementation of \kbd{rnfkummer},
these degrees must be \emph{prime}.) What about other extensions of degree
$2$ for instance?
\bprog
  L0= subgrouplist(bnr, [2])
  L = subgrouplist(bnr, [2], 1)
@eprog\noindent
\kbd{L0}, resp.~\kbd{L} is the list of those subgroups of the full ray class
group mod $7$, whose index is $2$, and whose conductor is $7$,
resp.~arbitrary. (Subgroups are given by a matrix of generators, in terms of
\kbd{bnr.gen}.) \kbd{L0} has $2$ elements, associated to the $2$ extensions
we already know. \kbd{L} has $7$ elements, the $2$ from \kbd{L0}, and
$5$ new ones:
\bprog
  L1 = eval( setminus(Set(L), Set(L0)) )
@eprog\noindent
The conductors are
\bprog
  vector(#L1, i, bnrconductor(bnr, L1[i]))
@eprog\noindent
among which one sees the identity matrix, i.e. the trivial ideal. (It is
\kbd{L1[3]} in my session, maybe not in yours. Take the right one!) Indeed,
the class group was cyclic of order $4$ and there exists a unique unramified
quadratic extension. We could find it directly by recomputing a \kbd{bnr}
with trivial conductor, but we can also use
\bprog
  rnfkummer(bnr, L1[3]) \\ @com pick the subgroup with trivial conductor!
@eprog\noindent
directly which outputs the (unique by Takagi's theorem) class field
associated to the subgroup \kbd{L1[3]}. In fact, it is of the form
$K(\sqrt{-\epsilon})$. We can check this directly:
\bprog
  rnfconductor(bnf, x^2 + bnf.fu[1])
@eprog\noindent

\subsec{Galois theory over $\Q$}

PARI includes a nearly complete set of routines to compute with Galois
extensions of $\Q$. We start with a very simple example.
Let $\zeta$ a $8$th-root of unity and $K=\Q(\zeta)$. The minimal
polynomial of $\zeta$ is the 8$th$ cyclotomic polynomial, namely
\kbd{polcyclo(8)} (=$x^4+1$).

We issue the command
\bprog
G = galoisinit(x^4 + 1);
@eprog\noindent
to compute $G=\text{Gal}(K/\Q)$.  The command \kbd{galoisisabelian(G)}
returns \kbd{[2,0;0,2]} so $G$ is an abelian group, isomorphic to $(\Z/2\Z)^2$, generated by
$\sigma$=\kbd{G.gen[1]} and $\tau$=\kbd{G.gen[2]}. These automorphisms are
given by their actions on the roots of $x^4+1$ in a suitable $p$-adic
extension. To get the explicit action on $\zeta$, we use
\kbd{galoispermtopol(G,G.gen[i])} for $i=1,2$ and get $\sigma(\zeta)=-\zeta$
and $\tau(\zeta)=\zeta^3$. The last non-trivial automorphism is
$\sigma\tau$=\kbd{G.gen[1]*G.gen[2]} and we have
$\sigma\tau(\zeta)=-\zeta^3$. (At least in my version, yours may return a
different set of generators, rename accordingly.)

We compute the fixed field of $K$ by the subgroup generated by $\tau$ with
\bprog
galoisfixedfield(G, G.gen[2], 1)
@eprog\noindent
and get $x^2 + 2$. Now we want the factorisation of $x^4+1$ over that
subfield. Of course, we could use \kbd{nffactor}, but here we have a much
simpler option: \kbd{galoisfixedfield(G, G.gen[1], 2)} outputs
\bprog
[x^2 + 2, Mod(x^3 + x, x^4 + 1), [x^2 - y*x - 1, x^2 + y*x - 1]]
@eprog\noindent
which means that
$x^4+1=(x^2-\alpha\*x-1)(x^2+\alpha\*x-1)$ where $\alpha$ is a root of $x^2+2$,
and more precisely, $\alpha=\zeta^3+\zeta$. So we recover the well-known
factorisation:

$$x^4+1=(x^2-\sqrt{-2}\*x-1)(x^2+\sqrt{-2}\*x-1)$$

For our second example, let us take the field $K$ defined by the polynomial
\bprog
  P = x^18 - 3*x^15 + 115*x^12 + 104*x^9 + 511*x^6 + 196*x^3 + 343;
  G = galoisinit(P);
@eprog\noindent Since \kbd{galoisinit} succeeds,
the extension $K/\Q$ is Galois. This time \kbd{galoisisabelian(G)} return
$0$, so the extension is not abelian, however we can still put a name on the
underlying abstract group. Use \kbd{galoisidentify(G)}, which return $[18,
3]$. By looking at the GAP4 classification we find that $[18, 3]$ is
$S_3\times\Z/3\Z$. This time, the subgroups of $G$ are not obvious,
fortunately we can ask PARI : \kbd{galoissubgroups(G)}.

Let us look for a polynomial $Q$ with the property that $K$ is the splitting
field of $Q(x^2)$. For that purpose, let us take $\sigma$=\kbd{G.gen[3]}.  We
check that \kbd{G.gen[3]\^{}2} is the identity, so $\sigma$ is of order $2$. We now compute the fixed field $K^\sigma$ and the relative factorisation of $P$ over
$K^\sigma$:
\bprog
F = galoisfixedfield(G, G.gen[3], 2);
@eprog\noindent
So $K$ is a quadratic extension of $K^\sigma$ defined by the polynomial
\kbd{R=F[3][1]}. It is well-known that $K$ is also defined by $x^2-D$
where $D$ is the discriminant of $R$ (over $K^\sigma$).
To compute $D$ we issue:
\bprog
D = poldisc(F[3][1]) * Mod(1,subst(F[1],x,y));
@eprog\noindent
Note that since \kbd{y} in \kbd{F[3][1]} denotes a root of \kbd{F[1]}, we
have to use \kbd{subst(,x,y)}.  Now we hope that $D$ generate $K^\sigma$ and
compute \kbd{Q=charpoly(D)}. We check that $Q=x^9+270\*x^6+12393\*x^3+19683$ is
irreducible with \kbd{polisirreducible(Q)}. (Were it not the case, we would
multiply $D$ by a random square.) So $D$ is a generator of $K^\sigma$ and
$\sqrt{D}$ is a generator of $K$. The result is that $K$ is the splitting
field of $Q(x^2)$.  We can check that with
\kbd{nfisisom(P,subst(Q,x,x\^{}2))}.

\section{Plotting}

PARI supports high and low-level graphing functions, on a variety of output
devices: a special purpose window under standard graphical environments (the
\kbd{X Windows} system, Mac OS X, Microsoft Windows), or a \kbd{PostScript}
file ready for the printer. These functions use a multitude of flags, which
are mostly power-of-2. To simplify understanding we first give these flags
symbolic names.
\bprog
  /* Relative positioning of graphic objects: */
  nw       = 0;  se       = 4;  relative = 1;
  sw       = 2;  ne       = 6;

  /* String positioning: */
  /* V */ bottom  =  0;   /* H */  left   = 0;   /* Fine tuning */ hgap = 16;
          vcenter =  4;            center = 1;                     vgap = 32;
          top     =  8;            right  = 2;
@eprog\noindent
We also decrease drastically the default precision.
\bprog
  \p 9
@eprog\noindent
This is very important, since plotting involves calculation of functions at
a huge number of points, and a relative precision of 28 significant digits
is an obvious overkill: the output device resolution certainly won't reach
$10^{28} \times 10^{28}$ pixels!

Start with a simple plot:
\bprog
  ploth(X = -2, 2, sin(X^7))
@eprog\noindent
You can see the limitations of the ``straightforward'' mode of plotting:
while the first several cycles of \kbd{sin} reach $-1$ and $1$, the cycles
which are closer to the left and right border do not. This is understandable,
since PARI is calculating $\sin(X^7)$ at many (evenly spaced) points, but
these points have no direct relationship to the ``interesting'' points on
the graph of this function.  No value close enough to the maxima and minima
are calculated, which leads to wrong turning points in the graph. To fix
this, one may use variable steps which are smaller where the function varies
rapidly:
\bprog
  ploth(X = -2, 2, sin(X^7), "Recursive")
@eprog\noindent
The precision near the edges of the graph is much better now.
However, the recursive plotting (named so since PARI subdivides intervals
until the graph becomes almost straight) has its own pitfalls.  Try
\bprog
  ploth(X = -2, 2, sin(X*7), "Recursive")
@eprog\noindent The graph looks correct far away, but it has a straight
interval near the origin, and some sharp corners as well.  This happens
because the graph is symmetric with respect to the origin, thus the middle 3
points calculated during the initial subdivision of $[-2,2]$ are exactly on
the same line.  To PARI this indicates that no further subdivision is needed,
and it plots the graph on this subinterval as a straight line.

There are many ways to circumvent this.  Say, one can make the right limit
2.1.  Or one can ask PARI for an initial subdivision into 16 points instead
of default 15:
\bprog
  ploth(X = -2, 2, sin(X*7), "Recursive", 16)
@eprog\noindent
All these arrangements break the symmetry of the initial subdivision, thus
make the problem go away.  Eventually PARI will be able to better detect such
pathological cases, but currently some manual intervention may be required.

The function \kbd{ploth} has some additional enhancements which allow
graphing in situations when the calculation of the function takes a lot of
time.  Let us plot $\zeta({1\over 2} + it)$:
\bprog
  ploth(t = 100, 110, real(zeta(0.5+I*t)), /*empty*/, 1000)
@eprog\noindent
This can take quite some time.  (1000 is close to the default for many
plotting devices, we want to specify it explicitly so that the result does
not depend on the output device.)  Try the recursive plot:
\bprog
  ploth(t = 100, 110, real(zeta(0.5+I*t)), "Recursive")
@eprog\noindent
It takes approximately the same time.  Now try specifying fewer points,
but make PARI approximate the data by a smooth curve:
\bprog
  ploth(t = 100, 110, real(zeta(0.5+I*t)), "Splines", 100)
@eprog\noindent
This takes much less time, and the output is practically the same.  How to
compare these two outputs?  We will see it shortly.  Right now let us plot
both real and complex parts of $\zeta$ on the same graph:
\bprog
  f(t) = my(z = zeta(0.5+I*t)); [real(z),imag(z)];
  ploth(t = 100, 110, f(t), , 1000)
@eprog\noindent (Note the use of the temporary variable \kbd{z}; \kbd{my}
declares it local to the function's body.)

Note how one half of the roots of the real and imaginary parts coincide.
Why did we define a function \kbd{f(t)}?  To avoid calculation of
$\zeta({1\over2} + it)$ twice for the same value of t.  Similarly, we can
plot parametric graphs:
\bprog
  ploth(t = 100, 110, f(t), "Parametric", 1000)
@eprog\noindent In that case (parametric plot of the real and imaginary parts
of a complex function), we can also use directly
\bprog
  ploth(t = 100, 110, zeta(0.5+I*t), "Complex", 1000)
  ploth(t = 100, 110, zeta(0.5+I*t), "Complex|Splines", 100)
@eprog

If your plotting device supports it, you may ask PARI to show the points
in which it calculated your function:
\bprog
  ploth(t = 100, 110, f(t), "Parametric|Splines|Points_too", 100)
@eprog

As you can see, the points are very dense on the graph.  To see some crude
graph, one can even decrease the number of points to 30.  However, if you
decrease the number of points to 20, you can see that the approximation to
the graph now misses zero.  Using splines, one can create reasonable graphs
for larger values of t, say with
\bprog
  ploth(t = 10000, 10004, f(t), "Parametric|Splines|Points_too", 50)
@eprog

How can we compare two graphs of the same function plotted by different
methods?  Documentation shows that \kbd{ploth} does not provide any direct
method to do so.  However, it is possible, and even not very complicated.

The solution comes from the other direction.  PARI has a power mix of high
level plotting function with low level plotting functions, and these functions
can be combined together to obtain many different effects.  Return back
to the graph of $\sin(X^7)$.  Suppose we want to create an additional
rectangular frame around our graph.  No problem!

First, all low-level graphing work takes place in some virtual drawing
boards (numbered from 0 to 15), called ``rectangles'' (or ``rectwindows'').
So we create an empty ``rectangle'' and name it rectangle 2 (any
number between 0 and 15 would do):
\bprog
  plotinit(2)
  plotscale(2, 0,1, 0,1)
@eprog
This creates a rectwindow whose size exactly fits the size of the output
device, and makes the coordinate system inside it go from 0 to 1 (for both
$x$ and $y$). Create a rectangular frame along the boundary of this rectangle:
\bprog
  plotmove(2, 0,0)
  plotbox(2, 1,1)
@eprog
Suppose we want to draw the graph inside a subrectangle of this with upper
and left margins of $0.10$ (so 10\% of the full rectwindow width), and
lower and top margins of $0.02$, just to make it more interesting. That
makes it an $0.88 \times 0.88$ subrectangle; so we create another rectangle
(call it 3) of linear size 0.88 of the size of the initial rectangle and
graph the function in there:
\bprog
  plotinit(3, 0.88, 0.88, relative)
  plotrecth(3, X = -2, 2, sin(X^7), "Recursive")
@eprog
(nothing is output yet, these commands only fills the virtual drawing
boards with PARI graphic objects). Finally, output rectangles 2 and 3 on
the same plot, with the required offsets (counted from upper-left corner):
\bprog
  plotdraw([2, 0,0,  3, 0.1,0.02], relative)
@eprog
\noindent The output misses something comparing to the output of
\kbd{ploth}: there are no coordinates of the corners of the internal
rectangle.  If your output device supports mouse operations (only
\kbd{gnuplot} does), you can find coordinates of particular points of the
graph, but it is nice to have something printed on a hard copy too.

However, it is easy to put $x$- and $y$-limits on the graph.  In the
coordinate system of the rectangle 2 the corners are $(0.1,0.1)$,
$(0.1,0.98)$, $(0.98,0.1)$, $(0.98,0.98)$.  We can mark lower $x$-limit by
doing
\bprog
  plotmove(2, 0.1,0.1)
  plotstring(2, "-2.000", left+top+vgap)
@eprog\noindent
Computing the minimal and maximal $y$-coordinates might be trickier, since
in principle we do not know the range in advance (though for $\sin(X^7)$ it
is easy to guess!). Fortunately, \kbd{plotrecth} returns the $x$- and
$y$-limits.

Here is the complete program:
\bprog
  plotinit(3, 0.88, 0.88, relative)
  lims = plotrecth(3, X = -2, 2, sin(X^7), "Recursive")
  \p 3          \\ @com $3$ significant digits for the bounding box are enough
  plotinit(2);      plotscale(2, 0,1, 0,1)
  plotmove(2, 0,0); plotbox(2, 1,1)
  plotmove(2, 0.1,0.1);
  plotstring(2, lims[1], left+top+vgap)
  plotstring(2, lims[3], bottom+vgap+right+hgap)
  plotmove(2, 0.98,0.1); plotstring(2, lims[2], right+top+vgap)
  plotmove(2, 0.1,0.98); plotstring(2, lims[4], right+hgap+top)
  plotdraw([2, 0,0,  3, 0.1,0.02], relative)
@eprog

We started with a trivial requirement: have an additional frame around
the graph, and it took some effort to do so.  But at least it was possible,
and PARI did the hardest part: creating the actual graph.
Now do a different thing: plot together the ``exact'' graph of
$\zeta({1/2}+it)$ together with one obtained from splines approximation.
We can emit these graphs into two rectangles, say 0 and 1, then put these
two rectangles together on one plot.  Or we can emit these graphs into one
rectangle 0.

However, a problem arises: note how we
introduced a coordinate system in rectangle 2 of the above example, but we
did not introduce a coordinate system in rectangle 3.  Plotting a
graph into rectangle 3 automatically created a coordinate system
inside this rectangle (you could see this coordinate system in action
if your output device supports mouse operations).  If we use two different
methods of graphing, the bounding boxes of the graphs will not be exactly
the same, thus outputting the rectangles may be tricky.  Thus during
the second plotting we ask \kbd{plotrecth} to use the coordinate system of
the first plotting.  Let us add another plotting with fewer
points too:
\bprog
plotinit(0, 0.9,0.9, relative)
plotrecth(0, t=100,110, f(t), "Parametric",300)
plotrecth(0, t=100,110, f(t), "Parametric|Splines|Points_too|no_Rescale",30);
plotrecth(0, t=100,110, f(t), "Parametric|Splines|Points_too|no_Rescale",20);
plotdraw([0, 0.05,0.05], relative)
@eprog

This achieves what we wanted: we may compare different ways to plot a graph,
but the picture is confusing: which graph is what, and why there are multiple
boxes around the graph?  At least with some output devices one can control
how the output curves look like, so we can use this to distinguish different
graphs.  And the mystery of multiple boxes is also not that hard to solve:
they are bounding boxes for calculated points on each graph.  We can disable
output of bounding boxes with appropriate options for \kbd{plotrect}.
With these frills the script becomes:
\bprog
plotinit(0, 0.9,0.9, relative)
plotpointtype(-1, 0)                \\@com set color of graph points
plotpointsize(0, 0.4)               \\@com use tiny markers (if available)
plotrecth(0, t=100,110, f(t), "Parametric|no_Lines", 300)
plotpointsize(0, 1)                 \\@com normal-size markers
plotlinetype(-1, 1)                 \\@com set color of graph lines
plotpointtype(-1, 1)                \\@com set color of graph points
opts="Parametric|Splines|Points_too|no_Rescale|no_Frame|no_X_axis|no_Y_axis";
plotrecth(0, t=100,110,f(t), opts, 30);
plotlinetype(-1, 2)                 \\@com set color of graph lines
plotpointtype(-1, 2)                \\@com set color of graph points
plotrecth(0, t=100,110,f(t), opts, 20);
plotdraw([0, 0.05,0.05], relative)
@eprog

\noindent Plotting axes on the second and third graph would not hurt, but
is not needed either, so we omit them.  One can see that the discrepancies
between the exact graph and one based on 30 points exist, but are pretty
small.  On the other hand, decreasing the number of points to 20 makes
quite a noticeable difference.

Keep in mind that \kbd{plotlinetype}, \kbd{plotpointtype},
\kbd{plotpointsize} may do nothing on some terminals.

Additionally, one can ask PARI to output a plot into a PS file: just
use the command \kbd{psdraw} instead of \kbd{plotdraw} in the above examples
(or \kbd{psploth} instead of \kbd{ploth}).  See \kbd{psfile} argument
to \kbd{default} for how to change the output file for this operation.  Keep
in mind that the precision of PARI PS output will be the same as for the
current output device.

Now suppose we want to join many different small graphs into one picture.
We cannot use one rectangle for all the output as we did in the example
with $\zeta({1/2}+it)$, since the graphs should go into different places.
Rectangles are a scarce commodity in PARI, since only 16 of them are
user-accessible.  Does it mean that we cannot have more than 16 graphs on
one picture?  Thanks to an additional operation of PARI plotting engine,
there is no such restrictions.  This operation is \kbd{plotcopy}.

The following script puts 4 different graphs on one plot using 2 rectangles
only, \kbd{A} and \kbd{T}:
\bprog
  A = 2;   \\@com accumulator
  T = 3;   \\@com temporary target
  plotinit(A);         plotscale(A, 0, 1, 0, 1)

  plotinit(T, 0.42, 0.42, relative);
  plotrecth(T, x= -5, 5, sin(x), "Recursive")
  plotcopy(T, 2, 0.05, 0.05, relative + nw)

  plotmove(A, 0.05 + 0.42/2, 1 - 0.05/2)
  plotstring(A,"Graph", center + vcenter)

  plotinit(T, 0.42, 0.42, relative);
  plotrecth(T, x= -5, 5, [sin(x),cos(2*x)], 0)
  plotcopy(T, 2, 0.05, 0.05, relative + ne)

  plotmove(A, 1 - 0.05 - 0.42/2, 1 - 0.05/2)
  plotstring(A,"Multigraph", center + vcenter)

  plotinit(T, 0.42, 0.42, relative);
  plotrecth(T, x= -5, 5, [sin(3*x), cos(2*x)], "Parametric")
  plotcopy(T, 2, 0.05, 0.05, relative + sw)

  plotmove(A, 0.05 + 0.42/2, 0.5)
  plotstring(A,"Parametric", center + vcenter)

  plotinit(T, 0.42, 0.42, relative);
  plotrecth(T, x= -5, 5, [sin(x), cos(x), sin(3*x),cos(2*x)], "Parametric")
  plotcopy(T, 2, 0.05, 0.05, relative + se)

  plotmove(A, 1 - 0.05 - 0.42/2, 0.5)
  plotstring(A,"Multiparametric", center + vcenter)

  plotlinetype(A, 3)
  plotmove(A, 0, 0)
  plotbox(A, 1, 1)

  plotdraw([A, 0, 0])
  \\ psdraw([A, 0, 0], relative)          \\ @com if hard copy is needed
@eprog

The rectangle \kbd{A} plays the role of accumulator, rectangle \kbd{T} is
used as a target to \kbd{plotrecth} only.  Immediately after plotting into
rectangle \kbd{T} the contents is copied to accumulator.  Let us explain
numbers which appear in this example: we want to create 4 internal rectangles
with a gap 0.06 between them, and the outside gap 0.05 (of the size of the
plot).  This leads to the size 0.42 for each rectangle.  We also
put captions above each graph, centered in the middle of each gap.  There
is no need to have a special rectangle for captions: they go into the
accumulator too.

To simplify positioning of the rectangles, the above example uses relative
``geographic'' notation: the last argument of \kbd{plotcopy} specifies the
corner of the graph (say, northwest) one counts offset from. (Positive
offsets go into the rectangle.)

To demonstrate yet another useful plotting function, design a program to
plot Taylor polynomials for a $\sin x$ about 0.  For simplicity, make the
program good for any function, but assume that a function is odd, so only
odd-numbered Taylor series about 0 should be plotted.  Start with defining
some useful shortcuts
\bprog
  xlim = 13;  ordlim = 25;  f(x) = sin(x);
  default(seriesprecision,ordlim)
  farray(t) = vector((ordlim+1)/2, k, truncate( f(1.*t + O(t^(2*k+1)) )));
  FARRAY = farray('t);  \\@com\kbd{'t} to make sure \kbd{t} is a free variable
@eprog\noindent
\kbd{farray(x)} returns a vector of Taylor polynomials for $f(x)$, which we
store in \kbd{FARRAY}.  We want to plot $f(x)$ into a rectangle, then make
the rectangle which is 1.2 times higher, and plot Taylor polynomials into the
larger rectangle.  Assume that the larger rectangle takes 0.9 of the final
plot.

First of all, we need to measure the height of the smaller rectangle:
\bprog
  plotinit(3, 0.9, 0.9/1.2, 1);
  opts = "Recursive | no_X_axis|no_Y_axis|no_Frame";
  lims = plotrecth(3, x= -xlim, xlim, f(x), opts,16);
  h = lims[4] - lims[3];
@eprog\noindent
Next step is to create a larger rectangle, and plot the Taylor polynomials
into the larger rectangle:
\bprog
  plotinit(4, 0.9,0.9, relative);
  plotscale(4, lims[1], lims[2], lims[3] - h/10, lims[4] + h/10)
  plotrecth(4, x = -xlim, xlim, subst(FARRAY,t,x), "no_Rescale");
@eprog

Here comes the central command of this example:
\bprog
  plotclip(4)
@eprog\noindent
What does it do?  The command \kbd{plotrecth(\dots, "no\_Rescale")} scales the
graphs according to coordinate system in the rectangle, but it does not pay
any other attention to the size of the rectangle.  Since \kbd{xlim} is 13,
the Taylor polynomials take very large values in the interval
\kbd{-xlim...xlim}.  In particular, significant part of the graphs is going
to be \emph{outside} of the rectangle. \kbd{plotclip} removes the parts of
the plottings which fall off the rectangle boundary
\bprog
  plotinit(2)
  plotscale(2, 0.0, 1.0, 0.0, 1.0)
  plotmove(2,0.5,0.975)
  plotstring(2,"Multiple Taylor Approximations",center+vcenter)
  plotdraw([2, 0, 0,  3, 0.05, 0.05 + 0.9/12,  4, 0.05, 0.05], relative)
@eprog\noindent
These commands draw a caption, and combine 3 rectangles (one with the
caption, one with the graph of the function, and one with graph of Taylor
polynomials) together. The plots are not very beautiful with the default
colors. See \kbd{examples/taylor.gp} for a user function encapsulating the
above example, and a colormap generator.

This finishes our survey of PARI plotting functions, but let us add
some remarks.  First of all, for a typical output device the picture is
composed of small colored squares (pixels), as a very large checkerboard.
Each output rectangle is a disjoint union of such squares.  Each drop
of paint in the rectangle will color a whole square in it.  Since the
rectangle has a coordinate system, it is important to know how this
coordinate system is positioned with respect to the boundaries of these
squares.

The command \kbd{plotscale} describes a range of $x$ and $y$ in the
rectangle.  The limit values of $x$ and $y$ in the coordinate system are
coordinates \emph{of the centers} of corner squares.  In particular,
if ranges of $x$ and $y$ are $[0,1]$, then the segment which connects (0,0)
with (0,1) goes along the \emph{middle} of the left column of the rectangle.
In particular, if we made tiny errors in calculation of endpoints of this
segment, this will not change which squares the segment intersects, thus
the resulting picture will be the same.  (It is important to know such details
since many calculations are approximate.)

Another consideration is that all examples we did in this section were
using relative sizes and positions for the rectangles.  This is nice, since
different output devices will have very similar pictures, while we
did not need to care about particular resolution of the output device.
On the other hand,
using relative positions does not guarantee that the pictures will be
similar.  Why?  Even if two output devices have the same resolution,
the picture may be different.  The devices may use fonts of different
size, or may have a different ``unit of length''.

The information about the device in PARI is encoded in 6 numbers: resolution,
size of a character cell of the font, and unit of length, all separately
for horizontal and vertical direction.  These sizes are expressed as
numbers of pixels.  To inspect these numbers one may use the function
\kbd{plothsizes}.  The ``units of length'' are currently used to calculate
right and top gaps near graph rectangle of \kbd{ploth}, and gaps for
\kbd{plotstring}.  Left and bottom gaps near graph rectangle are calculate
using both units of length, and sizes of character boxes (so that there
is enough place to print limits of the graphs).

What does it show?  Using relative sizes during plotting produces
\var{approximately} the same plotting on different devices, but does not
ensure that the plottings ``look the same''.  Moreover, ``looking the
same'' is not a desirable target, ``looking tuned for the environment''
will be much better.  If you want to produce such fine-tuned plottings,
you need to abandon a relative-size model, and do your plottings in
pixel units.  To do this one removes flag \kbd{relative} from the above
examples, which will make size and offset arguments interpreted this way.
After querying sizes with \kbd{plothsizes} one can fine-tune sizes and
locations of subrectangles to the details of an arbitrary plotting
device.

To check how good your fine-tuning is, you may test your graphs with a
medium-resolution plotting (as many display output devices are), and
with a low-resolution plotting (say, with \kbd{plotterm("dumb")} of gnuplot).

\section{GP Programming}

Do we really need such a section after all we have learnt so far ? We now
know how to write scripts and feed them to \kbd{gp}, in particular how to
define functions. It's possible to define \emph{member} function as well, but
we trust you to find them in the manual. We have seen most control
statements: the missing ones (\kbd{while}, \kbd{break}, \kbd{next},
\kbd{return} and various \kbd{for} loops) should be straightforward. (You
won't forget to look them up in the manual, will you ?)

Output is done via variants of the familiar \kbd{print} (to screen),
\kbd{write} (to a file). Input via \kbd{read} (from file), \kbd{input}
(querying user), or \kbd{extern} (from an external auxiliary program).
Perhaps the most useful simple command we haven't seen yet is
\kbd{allocatemem} which increase the size of \kbd{gp}'s ``scratch space''. If
you regularly see \kbd{PARI stack overflows!} messages, think about this
one.

To customize \kbd{gp}, e.g.~increase the default stack space or load your
private script libraries on startup, look up \kbd{The preferences file}
section in the User's manual. You probably want to use \kbd{trap} on startup
as well: it enables so-called \emph{break loops} when an error occurs, where
you get a chance, e.g to inspect (and save!) values of variables before the
prompt returns and all computations so far are lost. When break loops
are enabled, \kbd{Control-C} is treated as a non-fatal error: it freezes a
computation and gets you a new prompt so that you may e.g., increase debugging
level, inspect or modify variables (in fact, run arbitrary commands), before
letting the program go on.

For clarity, it is advisable  to declare local variables in user functions
(and in fact, with the smallest possible scope), as we have done so far with
the keyword \kbd{my}. One is usually better off avoiding global variables
altogether, but you are the sole master on board at this point.

To reach grandwizard status, you may need to understand the all powerful
\kbd{install} function, which imports into \kbd{gp} an (almost) arbitrary
function from the PARI library (and elsewhere too!), or how to use the
\kbd{gp2c} compiler and its extended types. But both are beyond the scope of
the present document.

Have fun !
\bye